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We report the non-perturbative tuning of parameters — Kc, Kb, and Kcrit — that are related to the bare heavy- 
quark mass in the Fermilab action. This requires the computation of the masses of D./' and Bs*-* mesons 
comprised of a Fermilab heavy quark and a staggered light quark. Additionally, we report the hyperflne splittings 
for d/' and Bi*^ mesons as a cross-check of our simulation and analysis methods. We find a splitting of 
145 ± 15 MeV for the Ds system and 40 ± 9 MeV for the Bs system. These are in good agreement with the 
Particle Data Group average values of 143.9 ± 0.4 MeV and 46.1 ± 1.5 MeV, respectively. The calculations are 
carried out with the MILC 2+1 flavor gauge configurations at three lattice spacings owO.15,0.12 and 0.09 fm. 



I. INTRODUCTION 



Lattice QCD calculations play a critical role in the study of standard model physics and the search for new physics. For a 
set of lattice QCD calculations to be viable, several basic tasks are necessary. The bare gauge coupling must be eUminated 
in favor of an observable allowing the conversion from lattice to physical units; the bare masses in the lattice action must be 
tuned to correspond to physical quarks; and experimentally established quantities must be calculated in order to substantiate the 
method's accuracy and rehability. Once these tasks are complete, a variety of quantities inaccessible to or not yet determined by 
experiment may be calculated, such as decay constants, form factors, and mass spectra. 

The Fermilab Lattice and MILC Collaborations have reported several calculations [1-8] based on ensembles of lattice gauge 
fields with 2+1 flavors of sea quarks, generated by the MILC Collaboration [9, 10]. Details of the scale setting can be found in 
Refs. [11, 12], and details of the Ught-quark mass tuning in Ref. [11]. In this paper, we report on the necessary tuning of the 
heavy-quark action for charmed and bottom quarks. In particular, we describe calculations of the heavy-light pseudoscalar and 
vector meson masses using, for light quarks, the asqtad staggered action [13] and, for heavy quarks, the Fermilab interpreta- 
tion [14] of the Sheikholeslami-Wohlert ("clover") action [15] for Wilson fermions [16]. We use the spin-average of these meson 
masses to nonperturbatively tune the hopping parameter k, which is equivalent to the bare heavy-quark mass. We also describe 
the determination of Kcrit, the value of k for which a degenerate Wilson pseudoscalar' s mass vanishes. The value of Kcrit plays a 
minor role in the calculation of heavy-light matrix elements [3-5], and a more important role when determining a renormalized 
quark mass [6]. Finally, as a by-product of these calculations, we report the spin-dependent hyperflne splittings for Bg and Dg 
mesons, which test how well we have improved the chromomagnetic interaction. 

Two aspects of the Fermilab method are important here. First, the Fermilab interpretation makes no assumptions about the 
size of the quark mass. Therefore, we are able to treat both charm and bottom quarks within the same framework. Second, 
since the Sheikoleslami-Wohlert action maintains the spin and flavor syimnetries of heavy quarks, heavy-quark effective theory 
(HQET) can be used to interpret and improve lattice discretization effects [17, 18]. HQET techniques can be used to show 
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how the improvement works for observables, such as meson masses, in a way simpler than, though equivalent to, the Symanzik 

improvement program [19]. 

This paper is organized as follows. Section II reviews the theoretical framework upon which these calculations are based. 
Section III contains specific descriptions of the gauge configurations, actions, and operators used for the meson masses. Sec- 
tion IV covers the components of the numerical analysis. Section V details the fitting procedures. Section VI presents the 
results for the non-perturbative tuning of the heavy-quark hopping parameters Kc and Kb, the hyperfine splitting, and the critical 
hopping parameter Kcrit- Section VII sunraiarizes with a discussion of improvements to these calculations that are currently 
underway. Details of the meson-mass discretization error estimation are given in Appendix A. Appendices B and C tabulate 
intermediate numerical results. The partially quenched chiral perturbation theory expression for the hyperfine spUtting is derived 
in Appendix D. 

II. THEORETICAL BACKGROUND 

The hopping-parameter form of the heavy-quark action is [14] 

S = So + Sb + Se, (2.1) 

where 

So = X]^"^" ~'^Y1 [^"(^ ~ 7^*)C^n,^V'n+A + ^n+/l(l +7/i)C^l,^V'n] , (2-2) 
n n,fj. 



Sb 



CbK ^ eijki^n^ijBn-k'tpn, 



(2.3) 



Se = icEK'^ilj^aoiEn^iipn, (2.4) 

where a^i, = 5 [7/^,7;/]. The chromomagnetic and chromoelectric fields Bn-i and En-i are standard and given in Ref. [14]. 
The term So includes dimension-five terms to alleviate the fermion doubling problem [16]. The couplings ce and cb of the 
dimension-five operators in Sb and Se are chosen to reduce discretization effects [14, 15]. 
The hopping parameter k is related to the tadpole-improved bare quark mass by 

amo = — (^--^], (2.5) 

where a is the lattice spacing, uq is the tadpole-improvement factor [20], and Kcrit is the value of k for which the pseudoscalar 
meson mass (of two degenerate Wilson quarks) vanishes. Our nonperturbative determination of Kcrit is discussed in Sec. VIC. 
To motivate our method of tuning k, we first discuss the meson dispersion relation. We then turn to the HQET description of our 
Lagrangian to understand how to best use the dispersion relation. 

The meson dispersion relation can be written, for |p| <C mo,a~^, as [14] 

E{p) = Ml + ^ + 0{p'). (2.6) 

Here, and throughout this work, we use lower-case m for quark masses and upper-case M for meson masses. Mi and M2 
are known as the rest mass and kinetic mass, respectively. Because the lattice breaks Lorentz invariance. Mi 7^ M2, although 
Ml — > M2 as a ^ for the action in Eq. (2.1). By tuning k, one could adjust the bare, heavy-quark mass such that either Mi 
or M2 is equal to the physical meson mass. (To set Mi = M2 requires the introduction, and tuning, of an additional parameter 
in the action. This is possible but, as discussed below, not necessary [14].) 

To clarify the role of the different masses in Eq. (2.6), it is useful to introduce an effective Lagrangian. This also sets up 
a language for discussing discretization errors later. Because the action in Eq. (2.1) has the same heavy-quark spin and flavor 
symmetries as continuum QCD, HQET is an obvious candidate for its description [17, 18]. To employ HQET, one separates the 
short-distance physics at the scale of the inverse heavy-quark mass 1 / mg from the long-distance physics at the characteristic 
scale of QCD, Aqcd- The fact that we have a lattice does not change the vahdity or utility of this separation. It simply means 
that the lattice spacing a must be included in the description of the short-distance physics. Thus, the short-distance coefficients 
of HQET applied to Eq. (2.1) differ from those arrived at by applying HQET to continuum QCD; these differences are the 
heavy-quark discretization errors. Parameters in the lattice action can be chosen to minimize them. 
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We introduce the heavy-quark effective Lagrangian for our lattice gauge theory by writing [17, 18] 

£lgt = Aight + £hqet, (2.7) 

where >Ciight is the Symanzik local effective Lagrangian for the light degrees of freedom and = means the Lagrangian on the 
right-hand side describes the on-shell matrix elements of the Lagrangian on the left-hand side. The HQET Lagrangian has a 
power-counting scheme, denoted by 

£hQET = ^^HQET' (2-8) 

s 

(s) 

where 'Chqet includes all operators of dimension 4 + s, with coefficients of dimension —s consisting of powers of the short 
distances, l/mg or a. The first few terms in £hqet are [17] 

4qet = -h^+\D4 + mi)h^+\ (2.9) 

£W = -h^+)^h^+) + -h(+)'^h^+), (2.10) 
^^^^ 2m2 2mB 

^'^^^ 8m| 8m% 

where is a two-component heavy-quark field, <r are the PauU matrices, and B and E are the continuum gauge fields. The 

masses mi,m2,mB,mE, and m/5 are functions of the bare-quark mass toq and the gauge coupling. For example, the masses 
nil and TO2 are defined to all orders in perturbation theory by Eq. (2.6), applied now to the pole energy of a one-quark state [21]. 

The entries in Eqs. (2.9)-(2.11) are commonly referred to as follows. £hqet gives the rest mass. The first term of £hqet i^ 

(2) 

the kinetic energy and the second is the chromomagnetic, or hyperfine, interaction. The first term of £hqet i^ spin-orbit 
interaction while the second is known as the Darwin term. 

For the pseudoscalar and vector meson rest masses, the HQET formalism can be used to show that [17] 

m[*^ =rm + Aut-^-dj^+ 0{l/m% (2.12) 
2m2 2mB 

where J is the total meson angular momentum with do = 3 and c?i = — 1 for the pseudoscalar (Mi) and vector (Mf ) mesons, 
respectively. The quantities Aiat, ^i.ut, and A2.iat are HQET matrix elements. At non-zero lattice spacing they contain dis- 
cretization effects from £iight, hence the subscript "lat". The continuum limit of these quantities yields their counterparts in 
HQET applied to continuum QCD [17], which provides a basis for computing the continuum-QCD quantities A and Ai [22]. 

Mass splittings and matrix elements such as decay constants and form factors are not affected by the value of mi [17]. Thus, 
Eqs. (2.9) and (2.10) show that the kinetic mass m2 is the first mass in the expansion that does play a role in the dynamics. We 
therefore would like to associate m2, and hence M2, with the physical mass, tolerating toi ^ m2 (and Mi 7^ M2) for nonzero 
lattice spacings. The nonperturbative tuning of k then entails adjusting k until the meson kinetic mass — determined by fits of 
Monte Carlo lattice data to the dispersion relation, Eq. (2.6) — equals that of the physical meson mass. A relation similar to 
Eq. (2.12) holds for M2 

M^*^ =1712+ Alat + 0(l/m), (2.13) 

with the leading discretization errors appearing in the 1/m contribution. Final values for the nonperturbative tuning of k are 
given in Sec. VI A. 
To calculate the hyperfine splitting of the Dg or Bg meson, consider 



FromEq. (2.12), 



Ai = Mi*-Mi. (2.14) 



Mi*-Mi=4^^-h---, (2.15) 
2m B 



which differs from the continuum spHtting only by discretization errors in the light quarks and gluons appearing in A2,iat. the 
mismatch of ms and its continuum counterpart (or, equivalently, the choice of cb), and similar contributions from higher- 
dimension operators [17, 23]. The splitting of kinetic masses, A2 = M| — M2, does not depend on m^; rather, it depends on 
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other generalized masses which are not tuned in our simulations.' Thus, Ai formally has smaller discretization errors than A2. 
Ai is also statistically cleaner than A2. In Eq. (2.15), l/ms is sensitive to the clover coupling cb in Eq. (2.3), so Ai tests how 
well it has been chosen. The Bg and Dg hyperfine splittings are given in Sec. VI B. 

III. SIMULATIONS 

In this section, we describe the gauge configurations used and the details of the actions, operators, and correlation functions 
that describe the heavy-light mesons. In Section HI A, we discuss the gauge configurations and the parameters that describe 
each ensemble. We also review how the lattice spacing is determined and the values of the conversion factors ri and r\/a. In 
Section in B, we discuss parameter choices for the valence quarks and the smearing of the heavy-quark wave function and how 
correlators are built from heavy and light quark fields. 

A. Gauge Configurations and Related Parameters 

We use the MILC gauge configurations [9, 10] that have 2+1 flavors of asqtad-improved staggered sea quarks [13] and a 
Symanzik-improved gluon action [24, 25]. Discretization errors from the sea quarks and gluons start at 0{asa?' , a^). The four- 
fold degeneracy of staggered sea quarks is removed by taking the fourth root of the determinant. To support the legitimacy of 
this procedure, Shamir has developed a renormalization-group framework for lattice QCD with staggered fermions, which he 
uses to argue that non-local effects of the rooted staggered theory are absent in the continuum Umit [26]. Additional support for 
this procedure comes from chiral perturbation theory arguments [27, 28]. Reviews of these papers and of other evidence that 
this procedure reproduces the correct continuum limit appear in [11, 29, 30]. 

Table I lists the parameters of the gauge configurations used in this work. All configurations have been gauge- fixed to Coulomb 
gauge. Ensembles of configurations are grouped by their approximate lattice spacing and are referred to as "fine" (a w 0.09 fm), 
"coarse" (a w 0.12 fm), and "medium-coarse" (a w 0.15 fm). The simulation bare masses of the light and strange sea quarks are 
denoted by amj and am'^, respectively, where amj is the mass of the two fighter sea-quarks. The range of amj is fight enough 
that the physical up- and down-quark masses can be reached by a chiral extrapolation, while am'^ is close to the physical strange- 
quark mass. For convenience below, we write (amj, am!^) to identify ensembles, e.g., "the (0.0031, 0.031) fine ensemble". Also 
in Table 1 are the tadpole factors uq [20, 31], determined from the mean plaquette and used to improve the gauge-configuration 
actions [9, 10]. The value of the physical strange-quark mass is denoted by the unprimed [31]. 

To convert between lattice and physical units, the physical value of the lattice spacing must be determined. We define the 
distance ri [12] by 

r?F(ri) = 1, (3.1) 

where F(r) is the force between static quarks, calculated on the lattice. For each ensemble, this yields a value of ri in lattice 
units, ri/a. The values are then "smoothed" by fitting ln(ri/a), from all ensembles, to a polynomial in /? and 2am'i+am'g [31]. 
The physical value of ri is obtained via the lattice calculation of an experimentally measurable quantity. We consider two current 
determinations here. One uses a lattice calculation of the T{2S)-T{1S) splitting [33] to arrive at ri = 0.318(7) fm [10, 34]. 
A more recent determination using ri/^ gives r\ = 0.3108(15)(t7g) fm [35]. These two determinations are consistent within 
errors. Because the determination of ri from f„ uses finer lattice spacings, we take that value, 

ri = 0.3108(t|g) fm (3.2) 

with no additional error. While this work was being completed, a new determination of ri that uses two mass splittings and one 

decay constant became available; ri = 0.3133(^3"^) [36], which is consistent with the value used in this work. Quantities can 
now be converted from lattice to physical units by using n and the appropriate value of ri /a given in Table I [31]. 

B. Meson Correlation Functions 

Table II lists the values of parameters used in the valence-quark actions. For the light valence quark, we again use the asqtad 
action [13] and masses am'g close to the physical value of the strange-quark mass, cf. Table I. From Eqs. (2.9)-(2.11), one can 



Tree-level expressions for these masses, and hence their mismatch, can be found in Ref. [23]. 
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TABLE I: Parameters describing the ensembles used. The dimensions of the lattice are given in terms of the spatial (Nl) and temporal (iVr) 
size in lattice units. The gauge coupUng is given by P = W/g^. The bare masses of the Ught and strange sea quarks are given by amj and 
am's , respectively. L — uNl is the linear spatial dimension of the lattice in fm. The column labeled A'^cf is the number of configurations 
used in this worlc. The plaquette-determined tadpole-improvement factor is uo [31]. The physical strange quark mass is arns [31] with errors, 
statistical and systematic, of less than one percent. The ratio ri /a is described in the text; errors are Hessian from the smoothing fit. The final 
column lists the value of the inverse lattice spacing using ri = 0.3108(lgg) fm to convert from ri/a; errors are from the error on ri and 
ri/a. 





N't X iVr 


13 


am'i 




L(fm) 


iVrf 


Wo 


arris 


n/a 


a-' GeV 


"Fine" a w 0.09 fm 


40^ X 96 


7.08 


0.0031 


0.031 


3.5 


435 


0.8779 


0.0252 


3.692(6) 


2.3441^^ 




28^ X 96 


7.09 


0.0062 


0.031 


2.4 


557 


0.8782 


0.0252 


3.701(5) 


2.349l!^J 




28^ X 96 


7.11 


0.0124 


0.031 


2.4 


518 


0.8788 


0.0252 


3.721(5) 


2.3621!;^ 


"Coarse" a w 0.12 fm 


24^ X 64 


6.76 


0.005 


0.05 


2.9 


529 


0.8678 


0.0344 


2.645(3) 


1.6791*^ 




20^ X 64 


6.76 


0.007 


0.05 


2.4 


836 


0.8678 


0.0344 


2.635(3) 


1.672l^^ 




20^ X 64 


6.76 


0.010 


0.05 


2.4 


592 


0.8677 


0.0344 


2.619(3) 


1.663lie 




20^ X 64 


6.79 


0.020 


0.05 


2.4 


460 


0.8688 


0.0344 


2.651(3) 


1.6831*^ 




20^ X 64 


6.81 


0.030 


0.05 


2.4 


549 


0.8696 


0.0344 


2.657(4) 


1.687li6 


"Medium-coarse" o w 0.15 fm 


16^ X 48 


6.572 


0.0097 


0.0484 


2.4 


631 


0.8604 


0.0426 


2.140(4) 


1.358l?^ 




16^ X 48 


6.586 


0.0194 


0.0484 


2.4 


631 


0.8609 


0.0426 


2.129(3) 


1.352tll 




16^ X 48 


6.600 


0.0290 


0.0484 


2.4 


440 


0.8614 


0.0426 


2.126(3) 


1.350tll 



see that with m2 timed to the physical mass, the leading mismatch between lattice and continutim physics is in the hyperfine term 
in >c[^Qgrp. In principle, one can tune tob to its continuum counterpart yielding a match between lattice and continuum actions 
for both terms in Eq. (2.10). Here, we use the tree-level expression for niB, which leaves the leading mismatch at 0{asaA). By 
setting ce = cb we obtain the Sheikholeslami-Wohlert, 0(a) -improvement of discretization errors in the action [15]. From the 
HQET perspective, this leaves ^ m2 in Eq. (2.11), but the effects of this mistuning are at O(a^A^) and 0{asaA? /rriQ). 
Implementing the improvements above and using tree- level tadpole improvement in the perturbative expressions [20, 24], we 
use Ce = Cb = ^• 

The values of uo used in the heavy-quark and light-valence actions are given in Table II. For the fine and medium-coarse 
ensembles, they are the plaquette values used to generate the MILC gauge configurations. For the coarse ensembles, the Landau- 
gauge link value was used. The use of different Uq definitions results in a slight mismatch between the light valence- and sea- 
quark actions. In part because the meson mass is relatively insensitive to the strange sea-quark mass, we do not expect any 
significant systematic errors from this mismatch. Changes in uq result in changes to the bare mass of the heavy quark as well, 
but this effect is partly absorbed by the nonperturbative tuning of n and Kcrit- Table II also lists the nominal values of the light 
valence-quark mass and sets of k values for bottom and charm mesons. These sets of k values, and mesons created from them, 
are referred to as charm-type or bottom-type. 

With the parameters of the actions set, we now turn to the construction of the two-point correlators. Contributions from excited 
states can be significantly reduced by using a spatially smeared source, sink, or both, for the heavy-quark propagator. For the 
correlators in this work, we use two types of source-sink combinations for the heavy quarks. One is simply a delta fimction for 
both the source and sink; we refer to this as the local correlator. The other smears the field "^{t, x) with a discretized version [37] 
of the 15 charmonium wavefunction, S{y), based on the Richardson potential [38]; 

ct>{t,x)=Y,S{y)ip{t,x + y), (3.3) 
y 

and the smearing wavefunction is applied after fixing to Coulomb gauge. Correlators using (^(t, x) are referred to as smeared 
correlators. All Ught valence quarks have a local source and sink. The meson correlator is 

CiAt^p) = Y.^O]{t,x) O,(0,0))e^P-, (3.4) 

where i.j denote the source, sink smearing of the heavy-quark field; for this work i = j. Oi{t, x) is a bilinear interpolating 
operator with a gamma-matrix structure that yields quantum numbers appropriate for either pseudoscalar or vector mesons. 
To construct this operator, we combine a one-component, staggered light-quark spinor with a four-component, Wilson-type 
heavy-quark spinor in a manner similar to Ref. [39], 

OE{t,x)=^^{t,x)r^^n^sit,x)xit,x), (3.5) 

where F = 75 or 7^; a, are spin indices; and n{x) = 7i ^72^73^74*. The fields tjj and x are the Wilson-type and staggered 
fields, respectively, and the smeared correlator is constructed in the same way, but with instead of tp. The transformation 
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TABLE II: Parameters used in the valence-quark actions. The bare masses of the light and strange sea quarks {am'i,am's) label the ensemble. 

The mass of the light (staggered) valence quark is given by am'^. ce and cb are the coefficients of the chromoelectric and chromomagnetic 
contributions to the Lagrangian. With ce ~ cb, they are the ususal Sheihkoleslami-Wohlert coupling, uo is the tadpole-improvement 
factor from measurements of the average plaquette for the fine and medium-coarse ensembles and from the Landau-gauge link on the coarse 
ensembles. Hopping parameter values n used for bottom-like and charm-like heavy quarks are given in the final two columns. 



Lattice 


{am'i, am's) 




Ce ~ Cb 




bottom-type k 


charm-type k 


Fine 


(0.0031, 0.031) 


0.0272, 0.031 


1.478 


0.8779 


0.0923 


0.127 




(0.0062, 0.031) 


0.0272, 0.031 


1.476 


0.8782 


0.090, 0.0923, 0.093 


0.1256, 0.127 




(0.0124, 0.031) 


0.0272, 0.031 


1.473 


0.8788 


0.0923 


0.127 


Coarse 


(0.005, 0.050) 


0.030, 0.0415 


1.72 


0.836 


0.086 


0.122 




(0.007, 0.050) 


0.030, 0.0415 


1.72 


0.836 


0.074, 0.086, 0.093 


0.119, 0.122, 0.124 




(0.010, 0.050) 


0.030, 0.0415 


1.72 


0.8346 


0.074, 0.086, 0.093 


0.119, 0.122, 0.124 




(0.020, 0.050) 


0.030, 0.0415 


1.72 


0.8369 


0.074, 0.086, 0.093 


0.122, 0.124 




(0.030, 0.050) 


0.030, 0.0415 


1.72 


0.8378 


0.086 


0.122 


Medium-coarse 


(0.0097, 0.0484) 


0.0387, 0.0484 


1.570 


0.8604 


0.070, 0.080 


0.115, 0.122^ 0.125 




(0.0194, 0.0484) 


0.0387, 0.0484 


1.567 


0.8609 


0.070, 0.076, 0.080 


0.115,0.122, 0.125 




(0.0290, 0.0484) 


0.0484 


1.565 


0.8614 


0.070, 0.080 


0.115, 0.125 



"Used only with am'^ = 0.484. 



properties of Ob. {x) under shifts by one lattice spacing are such that S can be viewed as playing the role of the (fermionic) taste 
index [30, 40]. In our correlation functions, Oh (a;) is summed over 2^ hypercubes, and so S can be interpreted as a taste degree 
of freedom in the sense of Refs. [41, 42]. 



IV. ANALYSIS OVERVIEW 



In this section, we describe the components of our analysis. Section IV A discusses the two-point correlator fits used to 
determine the meson energies aE(p). Section IV B describes how we fit the meson dispersion relation to obtain M2. Finally, 
Sec. IV C explains how k is tuned and how the hyperfine splitting is determined. 



A. Two-point Correlator Fits: E{p) 
To determine E{p), we simultaneously fit the local and smeared heavy-light-meson two-point correlators to the function 



Af-l 
ri=0 



Zf^ ^e-^^^P^* + e-EniPmT-t)'^ + (_l)t+l(^P^)2 ^g-E^ip)t ^ g-B^(p)(JVT-t) j 



(4.1) 



where Nt is the temporal extent of the lattice, and terms proportional to e~^'''^P)'^^^~*) are due to periodic boundary conditions. 
To simplify notation in this subsection, the lattice spacing a is not written out explicitly. Correlation functions containing 
staggered light quarks have contributions from both desired- and opposite-parity states with the opposite-parity states having the 
temporally-oscillating prefactor (—1)*+^ [39]. We take each energy level 77 in Eq. (4.1) to include a pair of states consisting of 
one desired- and one opposite-parity state; the number of pairs of states in a fit is given by N. Quantities associated with the 
tower of opposite-parity states are denoted by the superscript "p." 

Equation (4. 1) contains 2N exponentials, and the number of time slices in our data set is finite. Although it is straightforward 
to separate the two different parities — ^because of the (— 1)'+^ — it is difficult to separate states within each tower. Rather than 
relying solely on taking t large enough, we use the technique of constrained curve fitting [39, 43, 44]. We thus minimize an 
augmented [43], 

Xaug = X + / J ~2 ' (4.2) 

k Pk 

which means each fit parameter Pk is provided a prior Gaussian probability distribution function with central value and width 
{Pk, <^p^)- The central value for fitted quantities comes from minimizing x^^g on the whole ensemble. We take the parameters to 

be e'^^^ , ln(z|P^ ), and (for 77 > 0) ln(A£;^P' ), where AE^^^ = E^^^ - £;^p\, thereby enforcing a tower of states with increasing 
energy. 
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In general, one considers a quantity to be determined by the data only if the statistical error, discussed next, is smaller than the 
corresponding prior width. In this work, we are most concerned with the lowest-lying desired parity state, and the data — not the 
priors — always determine Eq and Z^ q. For parameters that are poorly constrained by the data, such as those describing excited 
states, these priors prevent the fitter from searching fruitlessly along flat directions in parameter space. Because of the freedom 
in choosing the prior, we test whether the ground-state results are prior-indpendent, and stable. When testing the stability of fit 
results, we use the Hessian error, defined as 




(4.3) 



because its straightforward definition allows it to be quickly calculated for a single fit. 

When using xtug, ^ measure the goodness of fit, we count the degrees of freedom as the number of data points; the number 
of fit parameters is not subtracted since there are an equal number of extra terms in xlug- ^ some cases, this could result in 
misleadingly low values of Xaug/'^^f- example, if the prior width ap^ is much larger than (Pk — Pk), the associated term in 
Xaug will be much smaller than the others. This could be adjusted a posteriori by reducing the degrees of freedom, but it would 
require devising a criterion for "large crp^". We do not make such adjustments in our analyses. Instead, to determine goodness 
of fit, we monitor the values of Xaug/'^of from constrained fits, but rely equally on the stability of fit results. 

We estimate statistical uncertainties by generating pseudo-ensembles via the bootstrap method. When fitting a pseudo- 
ensemble, the central value of each prior is drawn randomly from its Gaussian probability distribution while the prior width 
is kept the same [39, 43]. To prevent large, simultaneous but uncorrelated fluctuations among prior central values, which could 
destabilize a fit, we restrict the randomized prior central values to ±1.5(Tp. Final errors quoted for meson energies and functions 
thereof, such as the spin-averaged mass, are obtained from their bootstrap distributions. We define the upper (lower) 68%- 
distribution point as the value at which 16% of the distribution has a higher (lower) value. We refer to half of the distance 
between these two points as the average 68% bootstrap error. 



B. Dispersion Relation Fits: Tlie Kinetic Mass 

Having determined E(p), we use the dispersion relation to determine the kinetic meson mass, which we then use to tune the 
hopping parameter k. The low-momentum expansion for E{p) is 

i ^ 

where W4 and the deviation of M4 from M2 capture lattice artifacts. (In the continuum limit a^W^ = and M4 = M2.) The 
vector n is defined by 

ap = {2t, /Nl) n, (4.5) 

where is the spatial extent of the lattice, given in Table I; data are generated for \n\ < 3. Noise in E{p) increases with 
increasing momentum, though, and is substantial by the time 0{p'^) effects become significant. For charm-type mesons, squaring 
the energy yields a substantial cancellation in the ©(p^) contribution because aMi « aM2 ~ aM4. While this is not true for 
bottom-type mesons, the mass of these mesons is large enough to cause suppression via the 1/M factors whether E{p) or 
E'^{p) is used. By fitting to E'^{p) then, the contributions from O(p^) effects are reduced, and we are able to do a linear fit to 
low-momentum data, \n\ < 2. Setting Mi = £^(0) from the zero-momentum correlator, we square Eq. (4.4) and fit 

E^{p)-Ml =Cp'^ (4.6) 

to obtain C. Finally, we set M2 = Mi / C. The largest p is chosen so that the O(p^) effects are expected to be negligible, based 
on tree-level values of the analogous quark quantities W4 and l/m\. We confirm the negligibility of these terms by inspecting 
plots of the data and monitoring x^/dof. (We do not use constrained curve fitting here and so we minimize the usual x^-) This 
procedure is repeated for each bootstrap-generated pseudo-ensemble, yielding bootstrap distributions for aMi and aM2. 



C. The Hopping Parameter k and the Hyperflne Splitting Ai 



For tuning k, it is helpful to remove the leading discretization errors from spin-dependent terms. Let the spin-averaged kinetic 
meson mass be 

M2 = ^(M2+3M*), (4.7) 
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where M2 and M| are determined as described in Sec. IV B. This leaves the second, spin-independent term in Eq. (2.11) as the 
leading source of discretization error at O(a^A^). Our goal then is to determine the value of k that will result in a value of M2 
that agrees with the experimental value taken from the Particle Data Group (PDG). 

For each lattice spacing, we use the following procedure to tune k. Using three or more ensembles, we study the light sea- 
quark mass dependence of aM2 for at least one combination of k and m^. This gives us some insight into the behavior of 

aM2 in the physical-sea-quark-mass Umit and allows us to assign an uncertainty to 0M2 due to non-physical sea-quark masses. 
Next, on at least one ensemble, we determine a,M2 at two staggered, valence-quark masses near the strange-quark mass. This 
allows us to determine the dependence of aM 2 on the staggered, valence-quark mass and interpolate linearly to the physical 
value if no simulated mass is close enough to the tuned strange-quark mass. Having dealt with the staggered-valence and light 
sea-quark masses, we take aM2 at the physical, strange valence-quark mass at two values of k and interpolate Unearly in k to 
the spin-averaged value of the meson masses, given by the Particle Data Group (PDG) [45], converted to lattice units with a 
from Table I. Finally, we combine the uncertainties in the tuned value of k from statistical and discretization errors in the meson 
mass, staggered-valence mass mistuning, non-physical sea-quark masses, and errors from the lattice-spacing conversion of the 
PDG mass. 

To determine the hyperfine splitting, we start with the results for Mi = E{0). For each lattice spacing, we use values of 
aAi at, or linearly interpolated to, the tuned charm and bottom k values. We then consider uncertainties from statistics, the 
tuning of k and arus, non-physical sea-quark masses, and discretization. The value of aAi on the fine lattice is taken as our 
central value and results on the coarse and medium coarse lattices are used in the error analysis. In the final value, we also 
include an uncertainty due to the conversion to physical units. 

V. FITTING DETAILS FOR E{p) , Mi , M2 

In this section, we describe the details of our fitting procedure for the meson energy E{p) and the meson rest and kinetic 
masses. Mi and M2. Our objective here is to document thoroughly our fitting procedures, including values for the priors, and 
tests. Readers who are more interested in a summary can skip to Sec. V C. 

Section V A discusses the parameters used in our two-point correlator fits for E{p) (Sec. V A 1) and the evaluation of goodness 
of fit via Xaug/'iof tests of stability (Sec. V A 2). In most tests discussed here, Hessian errors were used, because they are fast 
and straightforward. Our complete data set, exhibited in Table II, contains several ensembles at each of the three lattice spacings. 
As explained in Sec. V A 1, one ensemble at each lattice spacing is chosen for the purpose of setting priors in Eq. (4.2). For 
tuning K, we need data over a range of k and am'^ on a fixed ensemble. At the fine lattice spacing, such data were generated on 
only one ensemble, (0.0062, 0.03 1), so we set priors and tune k on that same ensemble. For the coarse and medium-coarse lattice 
spacings, we have data for a range of k and am^ on several ensembles. We take the coarse (0.010, 0.050), and medium-coarse 
(0.0194, 0.0484) ensembles to set priors and then the ensembles with the smallest amj (and a range of k and am'q) to tune k. We 
compute the hyperfine splittings from the same ensembles on which k was tuned. These choices are summarized in Table III. 
Data from other ensembles listed in Table n are used to estimate uncertainties. 

Fits of the dispersion relation to determine M2 from E{p) are comparatively simple, and Sec. VB provides details that may 
be of interest. 



A. Two-point fits: E{p), Mi 

The number of gauge configurations in each ensemble is given in Table I. To improve statistics, we generate data at four time 
sources on each of the fine and coarse gauge configurations and at eight time sources for medium-coarse configurations. We 
also average the correlator points C{t) and C{Nt — t). In order to reduce the effect of correlations between data points from 
sequential configurations, we bin the data by groups of A^bin configurations. Because fits for this project were done in concert 
with other projects, iVbin = 4 was adopted. Comparisons of results using iVtin = 2,4, and 6 on the ensembles used here show no 
significant change in the fit-result error bars or the bootstrap distributions. To account for correlations in the two-point correlator 



TABLE in: Specific ensembles used in steps of the analyses. Setting priors is discussed in Sec. V A L Stability and goodness-of-fit tests done 
for E(p) results are described in Sec. V A2. K-tuning and hyperflne-splitting results are given in Sees. VIA and VI B, respectively. 



setting priors E{p) tests, tuning k, and the hyperfine splitting Ai 



Lattice 



Fine (0.0062,0.031) 
Coarse (0.010, 0.050) 

Medium-coarse (0.0194, 0.0484) 



(0.0062, 0.031) 
(0.007, 0.050) 
(0.0097, 0.0484) 
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data, the fitter uses the normalized, data-sample covariance matrix as an estimate of the correlation matrix. This matrix is remade 
for each bootstrap sample. 

/. Priors, time ranges, N 

We consider the setting of priors for the ground state parameters, excited-state amplitudes, and energy splittings separately. 
Ground-state (j] = 0) parameters are well-determined by the data; thus, the ground-state priors can, and should, be negligibly 
constraining. In contrast, energy splittings and excited state amplitudes are not well determined by the data, and the related 
priors are chosen such that they put reasonable bounds on the parameters. The next paragraphs describe how the priors are set. 
Note that the same set of priors is used for all ensembles at a given lattice spacing, for all momenta in the range \n\ = to 2, 
and for all k and am.[.^ of a given meson type, e.g., charm pseudoscalars. The priors used are tabulated in Tables IV-VI. 

We use information from a subset of our data, one ensemble per lattice spacing, to set the priors for the two-point-correlator 
fits. This is necessary because we do not have enough external knowledge to set them independently. The ensembles used to 
help set the priors are listed in Table III. Other ensembles are statistically independent of these ensembles and so the prior 
information can be viewed as external to fits on those ensembles. If possible, though, we do not want to exclude any data from 
our analysis, including the ensembles used in the setting of priors. For this reason, our procedure for setting priors keeps the 
amount of information we take from these ensembles to a minimum. Specifically, for a parameter P, we use averages over 
ranges of parameters, like the momentum, for the prior central value P and chose prior widths ap that are broad enough to cover 
the expected results for an entire subset of fits; e.g., the same priors are used for fits with \n\ = to 2. 

To set ground-state priors, we first fit to large-time data with A'^ = 1 in order to get a general idea of the ground-state 
parameter values. We then set A/^ > 1 and fit correlators at low and high momenta to ascertain the range of values the ground 
state parameters may take. We set prior central values for the ground-state energy of the desired- and opposite-parity states, 
aEo{p) and aEQ{p), at about the midpoint of the range seen in these fits. 

To understand our logic for setting the prior widths for aEo{p) and aEQ{p), recall that we use a Gaussian distribution for 
the prior P with a width a p. We set a^p and cr^pp large enough so that results across the entire momentum range used in the 
analysis should fall well within the ^-cr^p^, or 1-cr^^p, range of the distribution. After priors for the remaining parameters are 
set, we perform a complete set of fits and, for at least one ensemble at each lattice spacing, verify that, indeed, the final fit results 
for uEq and aE^ fit well within their respective prior distributions. 

Priors for the ground-state amplitudes are loosely based on the preUminary A'^ > 1 fits described above. In most cases, the 
central value is the nearest whole number to the average of these results. For the desired-parity state, the widths apaie chosen 
such that they easily span the range of values seen in the fits. For the opposite parity states, which are substantially noisier, the 
widths span the distance between the prior central value and the observed range in the results by about l-ap. 

Priors for all excited- state amplitudes were set to have a relatively small central value and a wide width. To set the prior for the 
energy splitting, we note that experimentally measured meson splittings are a few hundred MeV. We also bear in mind that the 
sum of a series of exponentials with a very small energy splitting is not a well-posed problem. Therefore, we chose the central 
value of the splitting to be several hundred MeV, slightly large, with a generous prior width. For example, on the fine lattice the 
prior for the splitting, ln(aA.E) = -1.45(1.0) is equivalent to AE w MeV. 

In the charm sector, the opposite-parity partner of the Ds{0^), the D*q{0^), is close to the DK threshold. In this case, the 
energy splitting should not be viewed as a meson mass splitting, and our choice of prior for the _D*q(0+) energy splitting may be 
inappropriate. The parity-partner signal is noisy, though, and in tests of the priors widths we see no change in the non-oscillating 
ground state energy aE{p), which is our main interest. For details, see Sec. V A2. 

To choose the time ranges for the fits, (imin. tmax), we first look at the data to determine the time by which the error in the 
data, e.g. the relative error in the correlator, has increased substantially. This gives us a potential value for tmax- From effective 
mass plots we can also see at what time slice the majority of the excited-state contamination has died off, giving us a potential 
value for imin- Constrained curve fitting is designed to reduce excited-state contamination of the lower-state fit parameters. 
Nevertheless, we do not see a significant reduction in the error from fitting to the smallest possible time slice, which requires 
including a larger number of states in the fit. For simplicity, we chose final time ranges that are the same for similar sets of data. 
These can be found in Table VII. 

With the time range set, we do fits for increasing values of the number of (pairs of) states N and look for the ground-state 
energy to stabilize. We choose the final values of N to be the minimum value needed to be in the stable region; these are given in 
Table VII. Figure 1 shows representative plots of aE{p) versus N from fits on the (0.0062, 0.031) fine ensemble. It is clear that 
for the minimum-value N, the central value of the fit result is always well within the stable region. In some cases, though, the 
(Hessian) error from the minimum- iV fit is smaller than that in the stable region. One could remedy this by choosing to fit with 
more states. Unfortunately, an increase in the number of states leads to non-gaussian bootstrap distributions with a significant 
number of outliers — clearly non-physical fit results that contain groimd states with low energies and very small amplitudes. 
Using the minimum possible number of states, no outliers have been seen in the distributions. 
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TABLE IV: Priors used for fine-ensemble two-point correlator fits for pseudoscalar and vector mesons. Priors for all higher amplitudes and 
splittings are the same as those for the first excited state. The fit-parameter numbers 15-20 label the second excited state and so on. A prior of 
ln(oAi5) = — 1.45l^ Q on the fine ensembles corresponds approximately to AE = 550I350 MeV. 



Charm Mesons Bottom Mesons 



fit parameter 


fit-parameter number pseudoscalar 


vector 


pseudoscalar 


vector 


Eo 


1 


0.90(40) 


0.90(40) 


1.75(60) 


1.75(60) 


K 


2 


1.0(40) 


0.95(40) 


1.85(60) 


1.85(60) 


ln(Zis,o) 


3 


1.0(2.0) 


1.0(2.0) 


1.0(3.0) 


1.0(3.0) 


4 


1.0(2.0) 


1.0(2.0) 


1.0(3.0) 


1.0(3.0) 


HZd,o) 


5 


-2.0(2.0) 


-2.0(2.0) 


-2.0(3.0) 


-2.0(3.0) 


6 


-2.0(2.0) 


-2.0(2.0) 


-2.0(3.0) 


-2.0(3.0) 


\n{AE) 


8 


-1.45(1.0) 


-1.45(1.0) 


-1.45(1.0) 


-1.45(1.0) 


In(ASP) 


9 


-1.45(1.0) 


-1.45(1.0) 


-1.45(1.0) 


-1.45(1.0) 


In(Zis.i) 
HZfs.i) 


10 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


11 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


ln(^d,i) 
HZlr) 


12 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


13 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 



TABLE V: Same as Table IV, but for the coarse ensembles. A prior of ln(aAi;) = -1.2lg:| on the coarse ensembles corresponds approxi- 
mately to AE = 500l^^^ MeV. 



Charm Mesons Bottom Mesons 



fit parameter 


fit-parameter number pseudoscalar 


vector 


pseudoscalar 


vector 


Eq 


1 


1.10(40) 


1.2(40) 


2.00(40) 


2.00(40) 


El 


2 


1.30(40) 


1.3(40) 


2.10(40) 


2.10(40) 


ln(^is,o) 
In(^fs.o) 


3 


1.0(2.0) 


1.0(2.0) 


1.0(2.0) 


1.0(2.0) 


4 


1.0(3.0) 


0.1(3.0) 


-1.0(2.0) 


-0.1(2.0) 


ln(^d,o) 

Mzio) 


5 


-1.0(2.0) 


-1.0(2.0) 


-2.0(2.0) 


-1.0(2.0) 


6 


-1.0(3.0) 


-2.0(3.0) 


-2.0(2.0) 


-2.0(2.0) 


\n{AE) 


8 


-1.2(0.5) 


-1.2(0.5) 


-1.2(0.5) 


-1.2(0.5) 


ln(A£;P) 


9 


-1.2(0.5) 


-1.2(0.5) 


-1.2(0.5) 


-1.2(0.5) 


ln(Zis,i) 
In(^rs.i) 


10 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


11 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


ln(^d,i) 

Hzi,,) 


12 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


13 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 



TABLE VI: Same as Table IV, but for the medium coarse ensembles. A prior of In(aAS) = -l.Otoj on the medium-coarse ensembles 
corresponds approximately to AE = SOOljoo MeV. 



Charm Mesons Bottom Mesons 



fit parameter 


fit-parameter number pseudoscalar 


vector 


pseudoscalar 


vector 


Eo 


1 


1.38(50) 


1.46(50) 


2.35(40) 


2.38(50) 


^0 


2 


1.50(60) 


1.58(60) 


2.48(50) 


2.50(50) 


HZisfi) 

HZfs.o) 


3 


0.48(1.0) 


0.95(1.0) 


0.12(1.4) 


0.60(1.0) 


4 


-0.65(1.0) 


0.20(1.0) 


-1.0(2.0) 


0.1(2.0) 


HZdfi) 
Hzio) 


5 


-0.90(1.0) 


-0.74(1.0) 


-1.15(1.0) 


-0.8(1.0) 


6 


-2.4(1.4) 


-1.8(2.0) 


-2.5(3.0) 


-1.8(3.0) 


ln{AE) 


8 


-1.0(0.5) 


-1.0(0.5) 


-1.0(0.5) 


-1.0(0.5) 


ln(A£;P) 


9 


-1.0(0.5) 


-1.0(0.5) 


-1.0(0.5) 


-1.0(0.5) 


In(^is.i) 
HZfs,i) 


10 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


11 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


In(Zd.i) 
In(^li) 


12 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


13 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 


-1.0(3.0) 
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FIG. 1 : Fitted values of aE{p) vs. the number of (pairs of) states for k = 0. 127, cliarm-type (a) pseudoscalar and (b) vector mesons and n = 
0.090, bottom-type (c) pseudoscalar and (d) vector mesons on the (0.0062, 0.031) fine ensemble. Results shown are for mesons with momenta 
n = (0, 0, 0) and (2, 0, 0). Errors are Hessian. 



2. Tests of Stability and Goodness-of-fit 

Having set the priors, time range, and number of states for the fits, we check the stability of the results and goodness of fit in 
several ways. For result stability, we check the effects of the time range used, the number of (pairs of) states N, and changes to 
the prior widths; we also compare the priors to the fit results. We look at a representative subset of fits for each lattice spacing: 
pseudoscalar and vector meson correlators at two different k values (one for charm and one for bottom) for a given light-valence 
mass, on one ensemble per lattice spacing, and with momenta n — (0, 0, 0) and (1, 1, 1) or (2, 0, 0). The specific values of k, 
am'q, and (amj, am'^) vary from test to test, and in some cases tests are extended to other values. A description of the data used 
in the tests discussed here can be found in Table Vlll. 



TABLE VII: Time range tmin-^max and number of (pairs of) states A'^ used in two-point correlator fits at each lattice spacing. For the time 
range, the first (second) number in parenthesis is fmin for the IS-smeared (local) correlator; imax is the same for both correlators. 



Lattice spacing 


Time range 


iV 


Fine 


(2,4)-25 


3 


Coarse 


(2,8)-15 


2 


Medium-coarse 


(5,6)-15 
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FIG. 2: Fit results shown as open (blue) circles are overlaid on the priors, black dots with dashed widths, for charm-type (a) pseudoscalar 
and (b) vector mesons and bottom-type (c) pseudoscalar and (d) vector mesons on the (0.0062, 0.031) fine ensemble, k = 0.127 and 0.090 
for charm- and bottom-type mesons, respectively; am'q = 0.0272. The upper [lower] plot is from a fit where the meson has momentum of 
n = (0, 0, 0) [(2, 0, 0)]. The fit-parameter numbers are defined in Table IV. In each panel, the leftmost cluster corresponds to quantities from 
the ground state; the middle cluster corresponds to the first excited state; and the right most cluster to the second excited state. Errors on the fit 
results are Hessian. For clarity, fit results are offset along the x-axis. 

For the time-range tests, we vary t^in over two to four time slices, increasing N if appropriate, and vary traax over five to ten 
time slices. We verify that there are no changes in the fit results beyond expected fluctuations.^ For number-of-states tests, we 



TABLE VIII: Data used in stability and goodness-of-fit tests. 



Lattice 



ensemble 



Fine (0.0062,0.031) 0.127; 0.090 or 0.093 0.0272 

Coarse (0.007,0.050) 0.122; 0.086 0.0415 

Medium-coarse (0.0097, 0.0484) 0.125; 0.070 0.0484 



^ In one case, k = 0.086, coarse (0.010, 0.005), although the ground-state energy is stable as tmax is varied, the value of x^/dof becomes large as tmax 
is increased beyond the final value (tmax = 15). This ensemble is not used directly for k tuning or hyperfine splitting determinations as explained in the 
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FIG. 3: Fit results shown as open (blue) circles are overlaid on the priors, black dots with dashed widths, for charm-type mesons on the (0.0062, 
0.031) fine ensemble, k = 0.127; am'q = 0.0272; n — (0, 0, 0). The upper plot is the same as the upper left (pseudoscalar) panel of Fig. 2 (a). 
The lower plot is from a fit which only differs by the use of = 4 pairs of states. The fit-parameter numbers are defined in Table IV. In each 
panel, the leftmost cluster corresponds to quantities from the ground state; the middle cluster corresponds to the first excited state; the next 
cluster corresponds to the second excited state and so on. The (desired-parity) ground-state quantities are stable to this change while other, 
excited-state, parameters are not. Errors on the fit results are Hessian. For clarity, fit results are offset along the i-axis. 



verify that the result is stable as N is increased. Figure 1 shows example results for the (0.0062, 0.031) fine ensemble. Similar 
results are seen for the coarse and medium-coarse ensembles and for the ground-state amplitudes Zig and Z^. 

For prior-width tests, we reduce the widths by a factor of two for the non-oscillating ground state quantities and the energy 
splittings and repeat the fits. All changes observed are within statistical errors and, in most cases, the changes are substantially 
smaller than one a. For charm, we also test for effects of the DK threshold near the £'*g(0+) state. This splitting is 50 to 100 
Me V, which is a several-cr^ ^ deviation from our prior central value. We ran separate tests on each lattice spacing using a prior 
width of ^ — 2.5 for the oscillating-state energy splitting. In units of MeV, this puts a 50-MeV splitting within la^ ^ 
of the prior central value. The ground and first-excited-state energies of the oscillating state are affected by this change but not 
in a systematic way. This indicates that the oscillating-state signal is not strong in our data. Our main interest, though, is the 
non-oscillating ground state energy aE{p); this value is unaffected by the change in cr^ ^ . 

In addition, we compare fit results with their priors. Figure 2 gives examples of these comparisons for fits on the (0.0062, 
0.031) fine ensemble for charm- and bottom-type mesons. The x-axis labels the fit-parameter number, defined in Table IV; the 
ground-state energy and amplitudes of the desired-parity state are at positions 1, 3, and 5. We find that fit results for ground-state 
quantities are well within the prior widths. For excited states, in some cases the fitter simply returns the prior value, indicating 
that the quantity is not constrained by the data. In other cases, the results appear to be constrained by the data, indicating that 
some excited-state signal is in the correlator and the fitter adjusts the amplitudes to absorb it. Although it may appear in Fig. 2 
that a number of excited-state quantities are well-determined, this is an artifact of a minimum-iV fit; unlike the ground-state 
parameters, the excited state results are not stable as N is increased. For example. Fig. 3 compares the fit results shown in 
the upper left (pseudoscalar) panel of Fig. 2 (a), which uses = 3, with a fit which only differs by the use of = 4. The 
comparison demonstrates that the (desired-parity) ground-state quantities are stable to the change in N while other, excited-state, 
parameters are not. 

For goodness-of-fit we begin by looking at the augmented x^/dof for each fit and verify that it is « 1 or smaller, where "w 1" 
is based on the 80% range of the x^/dof distribution for a given number of degrees of freedom. As a final check, we overlay the 
result on an effective-mass plot. We define the "effective energy" 

2aEcs{p) = In [C(t)/C(t + 2a)] (5.1) 

using a step of two time units in order to accommodate the oscillating contribution from the opposite-parity state. Figure 4 
shows plots comparing aEcs{p) to the fit result on the (0.0062, 0.031) fine ensemble. The ground-state-energy result from the 



introduction to this section. 
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FIG. 4: Effective energy plots, aEcsip), for charm-type (a) pseudoscalar and (b) vector mesons and bottom-type (c) pseudoscalar and (d) 
vector mesons on the (0.0062, 0.031) fine ensemble, k = 0.127 and 0.093 for charm- and bottom-type mesons, respectively; am'q = 0.0272. 
The upper [lower] plot is from a fit where the meson has momentum of n = (0, 0, 0) [n = (2, 0, 0)]. Open (blue) triangles mark the local 
correlator and open (red) circles mark the IS-smeared correlator. Lines connecting the data points are simply to guide the eye; they are not a 
fit. The unadorned black line is the multi-correlator fit result and the shaded band marks the average 68% bootstrap error. 



multiple-state fit is shown as a straight line segment over the time range fit. The band encompasses the average 68% bootstrap 
error. In each case, the fit result nicely matches the effective-energy plateau. 



B. The kinetic mass M2 

Given results for aE{p), we fit data where \n\ < \/3 to Eq. (4.6) to determine the pseudoscalar and vector kinetic meson 
masses. Fits use a correlation matrix constructed from the bootstrap distributions. The tables in Appendix B give results for aM2, 
aM|, and a A/2 on the ensembles used for tuning, listed in Table III. Included in the tables are the x^/dof and the probability 
that would exceed the value from the fit, known as the p value [45]. Typical dispersion relation fits are shown for the (0.0062, 
0.031) fine ensemble in Fig. 5. 

In addition to statistical eiTors, we consider uncertainties from unphysical sea-quark masses, mistuning of the valence strange 
quark, and discretization. The noise in M2 makes it difficult to discern how M2 depends on the sea-quark masses. The Mi data 
is much cleaner, though, and we can use it to estimate the sea-quark error on M2, and hence k. To do this, we first note that, 
cf. Eq. (2.12), 



aMi = ami + aAiat + 0(1/toq) 



(5.2) 
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FIG. 5 : Results of fits to the dispersion relation for (a) charm-type (k = 0. 1 27) and (b) bottom-type (k = 0.0923) mesons on the (0.0062, 0.03 1 ) 
fine ensemble. (Blue) dots are the data. A black line shows the fit result with the (pink) shaded band showing the one-sigma error from the fit. 
Upper panels show results for pseudoscalars and lower for vectors. 
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FIG. 6: The spin-averaged meson rest mass in physical units versus the ratio of the light to strange sea-quark masses mj / m's for the (a) fine and 
(b) coarse ensembles. Error bars are statistical only, from the average 68% bootstrap error. On the far right of the plot is a (red) bar indicating 
the size of the one-sigma statistical error on riM2- The upper-panel plot is for charm-type mesons, lower is for bottom-type. Values of k used 
are 0.127, 0.0923 on the fine ensembles with am'q = 0.0272 and k 0.122, 0.086 on the coarse ensembles with am'q = 0.0415. 



aMa = aTO2 +aAiat +0(l/mQ) (5.3) 

where ami and am2 capture the leading heavy-quark dependence and Ajat depends only on the light degrees of freedom. 
Taking aAiat to be the same for both aMi and aM2 (see Appendix A and Ref. [46]) we can estimate the size of the effect of 
non-physical (light) sea quark masses on aAiat, and hence aM2, by studying the behavior of aMi as the light sea-quark masses 
are varied. 

In Fig. 6, we plot the spin-averaged meson rest mass riM i versus the ratio of the light to strange sea-quark masses m[/m'g 
for the coarse and fine ensembles used here. On the far right of each plot is a bar indicating the size of the 1-ct statistical error on 
riM2\ for fine this is from the (0.0062, 0.031) ensemble and for coarse the (0.007, 0.050) ensemble. The light sea-quark mass 
dependence is negligible compared to the statistical error on riAfa- We find similar behavior for the medium-coarse ensemble. 

We must also consider how the non-physical value of the strange sea-quark mass affects M2. The strange sea-quark mass is 
mistuned by an amount G.lOam^, O.Slam^ and 0.12am^ on the fine, coarse, and medium-coarse ensembles, respectively. The 
continuum chiral perturbation theory expression for the heavy-light spin-averaged mass [47] shows that the leading sea-quark 
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dependence of M2 is proportional to the sum over the sea-quark masses, 2mJ + m',. Hence, varying amj tells us about the 
effect of varying am'^. Figure 6 shows that a change of O.SaWg in am'i has a neghgible effect on M2, so we conclude that the 
mistuning of am'^ has a negligible effect as well. 

The tuned value of the strange-quark mass on each ensemble is given in Table I. On the fine lattice, the valence-quark mass 
used in the simulation, am'^ = 0.0272, differs from the physical value anis = 0.0252 by 0.0020. A comparison of our results 

for aM2 in Table XVIII shows that even a deviation in om^ of twice this size does not discemibly affect aM2. The situation is 

similar for the coarse and medium-coarse results. For the coarse ensembles, the simulation mass am'^ = 0.03 differs by 0.0044 

from the tuned value of anig. Table XIX shows that aM2 is barely affected at the l-f^j^^^ level as am'g changes by over twice 
this size. For the medium-coarse ensembles, the simulation mass of 0.0484 differs from the tuned strange-quark mass by 0.0058. 
A comparison of the values of aM2 in Table XX shows that a deviation in am'^ just under twice this size yields, at most, a 

^'^a.M2 variation in aM2. Therefore, we take our results of aMa at am'g = 0.0272, 0.03, and 0.0484 as the masses of the Bs 
and Dg on the fine, coarse and medium-coarse ensembles, respectively, with no additional error for valence-mass mistuning. 

In Appendix A, we derive an expression for the discretization error in M2, M2 = Mcontinuum + 6M2. The result, Eq. (A22), 
can be written 



SMo 



6m2 



—^-1]+ 4w4{m2ay 



(5.4) 



replacing (p^) of Eq. (A22) with A^. Expressions for the short-distance coefficients m2, m4, and W4 are given in Ap- 
pendix A [14, 23]. To estimate the discretization error, we use values of the physical (pole) quark mass (1.4 GeV for charm and 
4.2 GeV for bottom) for m2 in the prefactor of Eq. (5.4), and A = 0.7 GeV. Using these values, uq from Table II, and Kcrit from 
Table XVII yields the values of SM2 shown in Table IX. The error estimate in Eq. (5.4) pertains to the kinetic mass, but the main 
focus here is the tuning of k. After tuning, we shall propagate this error from M2 to Kc and Kb- 



C. Fitting Summary 

The preceding subsections contain many details intended for those engaged in similar analyses. In this section, we re- 
emphasize the main features of the analysis. Because, in this and related [3-8] work, we are interested in the ground state, 
we do not dwell on the excited states here. 

Our priors are guided by the data, using one ensemble to set them and (generally) other ensembles for physical results. We 
choose a time range such that the fit results for the ground state are stable, listed in Table VII. We also test for stability as 
the number N of (pairs of) exponentials grows — as shown in one example in Fig. 1 — and choose the minimum value of N for 
which the central value is stable within errors. The errors on the ground-state amplitudes and energies are always determined 
by the data, not the priors, as shown in Fig. 2 and 3. (In many cases, even excited-state information is data-determined, not 
prior-determined.) Figure 4 shows that the fits agree with the effective energies. (Note that the oscillations of aEes at small t are 
to be expected with staggered quarks.) In conclusion, the constrained curve fitting for E{p) has worked as advertised, subsuming 
the subjectivity of fit ranges and different choices of N into robust results for both central value and error bar. Figures 5 and 6 
show that, once E{p) is well-determined, we can straightforwardly obtain the kinetic mass M2 and the hyperfine splitting. 



VI. RESULTS 



In this section, we present the main results of these calculations, including our error analysis. Section VIA focuses on the 
tuned values of Kc and Kb, Sec. VI B on the Dg and Bg hyperfine spUttings, and Sec. VI C on the critical value of the hopping 



TABLE IX: The relative error in the tuned hopping parameter 5k/ k, due to discretization effects in the kinetic meson mass. The ensembles 
used are (0.0062, 0.031), (0.007, 0.050), and (0.0097, 00484) for the fine, coarse, and medium-coarse lattices, respectively. Values of k are 
0.127 and 0.0923 on fine; 0.122 and 0.086 on coarse; and, 0.122 and 0.076 on medium-coarse. The [• • •] denotes the quantity in brackets in 
Eq. (5.4). We use (A^/6mch) = 0.0583 and (A^/6mbot) = 0.0194 to convert the [• • •] to SM2. Values of 6k/k are given as fractions not a 
percentage. 
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-0.0256 


coarse 
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parameter Kcrit- 



A. The tuning of Kc and Kb 

As discussed in Sec. VB, effects from non-physical sea-quark masses and the mistuning of the valence strange-quark mass 
are negligible compared to the statistical error on aM2- In that section, we explain why taking aM2 at certain values of am'^ is 
an acceptable approximation to aM2 at the tuned physical strange-quark mass. We choose to tune n at those same am'^, which 
are am'g = 0.0272 on the (0.0062, 0.031) fine ensemble, am'g = 0.03 on the (0.007, 0.050) coarse ensemble, and am'g = 0.0484 
on the (0.0097, 0.0484) medium-coarse ensemble. 

To obtain the tuned k for the charm (bottom) quark, Kc {k^), we want to interpolate M2 to the PDG value of the spin-averaged 
Dg (Bg) mass [45]. In practice, it is simpler to do the interpolation with the meson mass in lattice units. Hence, we linearly 
interpolate aM2 to aMpDc the PDG value for the meson mass converted to lattice units with a from Table I. This interpolation 
is repeated for the entire bootstrap distribution of aM2- We then estimate the statistical error on k as the average 68% bootstrap 
error described in Sec. IV A. The discretization error in M2, SM2, is given by Eq. (5.4), and is always positive. This results 
in a single-sided, negative error bar on k. We convert dM2 to the error, Sk, using dM2/dK ~ dm2/dK and expressions for 
moo and m20 given in Appendix A. The Sk are given in Table IX. The experimental errors on the PDG values are negligible. 
The remaining errors to consider are those which appear in the conversion between lattice and physical units. The error in the 
determination of ri/a is negligible, so we only need to consider the error in ri, given in Eq. (3.2). 

The error on n is propagated to an error on and then to an error on oMpdg> denoted ctpdg- Table X gives the values 
of the PDG meson masses used in this work and tabulates their spin-averaged mass and hyperfine splitting. Table XI gives the 
spin-averaged mass in lattice units. The uncertainty ctpdg is propagated to k using the standard error formula ctr = ctpdg / s, 
where s is the slope used in the interpolation. Table XII gives the error budget for and k^,, and Table XIII lists the final tuned 
results. 



TABLE X: PDG values of the pseudoscalar and vector masses for the Ds and Bs mesons and the hyperfine splitting A [45]. Also listed is the 
derived quantity M, the spin-averaged mass. 

M (GeV) M* (GeV) M (GeV) A (MeV) 
Ds 1.96849(34) 2.1f23(5) 2.0763(4) 143.9(4) 
Bs 5.3661(6) 5.4120(12) 5.4005(9) 46.1(1.5) 



TABLE XI: Spin-averaged PDG masses converted to lattice units with an error from the uncertainty in the lattice spacing a. Values of a used 
in the conversion can be found in Table I. 
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Medium-coarse (0.0097, 0.0484) 1.5291^:^3^ 3.9771^;' 



TABLE XII: Percent errors in the tuned k and the total error. For several sources of uncertainty, we determined that the error was smaller than 
the precision of these calculations. This is indicated by an entry of "0.0" in the table. 
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B. The rest mass and hyperflne splitting 



In this section, we discuss the uncertainties in our calculation of the hyperflne spUtting and compare our flnal results, for the 
Bs and Dg systems, with the PDG values. To support the discussion, we tabulate our results for the pseudoscalar and vector 
meson rest masses and the hyperflne splitting, aM, aM*,aAi, ri Ai, in Tables XXI-XXIII in Appendix C. Statistical errors 
in these tables are the average 68% bootstrap errors described in Sec. IV A. The other errors we consider are the mistuning of 
the valence strange-quark mass, unphysical sea-quark masses, the uncertainty in the tuning of k, discretization effects, and the 
conversion to physical units. For the central value, at each lattice spacing, we take aAi at the tuned values of and Kb, linearly 
interpolating in k when necessary. 

PDG results for the hyperflne splitting show a weak dependence on the Ught-quark valence mass, so we expect the mis- 
tuning in the simulated valence strange-quark mass to have a negligible effect.^ The simulation valence masses am^ = 
0.0272. 0.03, 0.0484 for the flne, coarse, and medium-coarse lattices, respectively, differ from the physical anig given in Ta- 
ble I by 0.0020, 0.0044, 0.0058, respectively. Tables XXI-XXIII show that, indeed, these small mistunings have a negligible 
effect on the hyperflne spUtting. Hence, we do not interpolate to anis; rather, we take aAi at the valence masses am'g Usted 
above as the result at the physical strange valence-quark mass and take the error for this approximation to be negligible. 

To estimate the error due to the non-physical values of the sea-quark masses we use partially-quenched chiral perturbation 
theory. The needed expression is derived in Appendix D and we repeat Eq. (Dl) here for convenience. The hyperflne splitting 
M* — Mx of a heavy-light meson with light- valence quark x is 



which must be determined from the lattice data. Working at a flxed value of rux, we can use the difference of spUttings at 
different values of m; to determine A^'^^ Given A''^^ we can flnd the difference between the spUtting at simulation values 
of {rn'i.m'g) and the physical values (m;,phys, mg^phys)- We take this difference as the error due to the non-physical sea-quark 
masses. 

We have tabulated values of the hyperflne splitting in physical units, ri Ai, in Appendix C 2. Figure 7 shows how ri Ai varies 
with the light sea-quark mass on flne and coarse lattices. From Fig. 7, it is clear that, due to statistical variation in the splitting, 
using the difference in the central values of splittings from any two points will yield different values for A*^"^). For the flne and 
coarse ensembles, we look only at the ami / o-'mg = 0.4 to 0. 1 and ami / amis = 0.4 to 0.2 differences and take the one that gives 
the larger error; for medium coarse, we have no ami/anis =0.1 data and so take the error from the ami /arris = 0-4 to 0.2 
difference. 

For the error estimate, we take / = 131 MeV and = 0.51 [48]. We relate meson to quark masses by 



where Bq is determined empirically with tiBq = 6.38, 6.23, 6.43 on the flne, coarse, and medium-coarse lattices, respectively. 
These values of Bq come from tree-level flts to MILC light-meson data, as described in Refs. [2, 11, 35]. We calculate A^*^) for 
each meson type, Bg and Dg, at each lattice spacing. We then calculate the difference 



where the subscript "sim" ("phys") denotes simulation (physical) sea-quark mass inputs {am,i,ams)- For the physical masses, we 
use (am;,phys,am5,phys) = (0.00092, 0.0252), (0.00125, 0.0344), (0.00154, 0.0426) for the flne, coarse, and medium-coarse 
lattices, respectively. These values of the quark masses are taken from Ref. [11], after adjustment for the ri scale used here. 




= Bo{mro+my) 



(6.2) 



(M* - M,),i^ - 



(M* - M,)phys 



(6.3) 



TABLE XIII: Final tuned results for Kc and Kb with the total error. 



Fine 



Coarse Medium-coarse 



Kc 0.127(2) 0.1219l|., 0.122i;^ 
Kb 0.090(5) 0.082(8) ().()77^-^ 



^ For X = B or D, the difference between the Mx* -Mx^ spUtting and the Mx»-Xx splitting is measured to be about 1% or less [45]. 
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The simulation masses are those on the the (0.0062, 0.031) fine, (0.007, 0.050) coarse, and (0.097, 0.0484) medium-coarse 
ensembles. The error calculated in this manner is labeled "sea-quark masses" in Tables XIV and XV. 

For the uncertainty in aAi due to the error in k, recall that the non-negligible sources of error in k, from Table XII in Sec. VI A, 
are statistics, units conversion, and discretization error in M2. Because we want to consider discretization errors separately from 
all others, we start by considering only the K-tuning error that comes from statistics and units-conversion. To convert the error 
in K to an error in aAi, we look at the change in aAibetween two values of k on the (0.0062, 0.031) fine, (0.007, 0.050) coarse, 
and (0.0097, 0.0484) medium-coarse ensembles; specific values can be found in Tables XXI-XXIII. This is the error labeled "k 
tuning" in Tables XIV and XV. 

For the Dg (Bg) meson. Table XIV (XV) gives the error budget for aAi at each lattice spacing, from all sources except 
discretization. These are statistics, valence-mass mistuning, unphysical sea-quark masses, and k tuning. In Fig. 8, these values 
are plotted as black, filled dots. 

We now consider the three, distinct sources of discretization error in aAi. The first is indirect, coming from the discretization 
error in aM2, which is propagated to an error on k as discussed in Sec. VIA. This error can be traced to a mismatch between 
the spin-independent O(p^) terms in Eq. (2.8) (not given explicitly) and the corresponding terms in the effective Lagrangian 
for continuum QCD. These terms contribute to aM2 as discussed in Appendix A. The second source of discretization error 
is a direct result of the lattice-continuum mismatch of the dimension-seven operator {icr ■ B, D^} [23].'^ The third source of 
discretization error is the 0{as) mismatch in the coefficient of the ia ■ B operator in Eq. (2.10). For the discussion of error 
estimates below, it is useful to recall that the heavy-quark dynamics associate m2 with the physical quark mass. Mismatches 
between m^. and the generalized masses associated with other operators capture the heavy-quark discretization effects. We now 
give numerical estimates of the error from each source. 

Our estimate of discretization error in aM^ and its inclusion in the error on k is discussed in Sec. VIA. In Fig. 8, the value of 
ri Ai with an error that includes only the uncertainty due to the discretization error on k is shown as an open (blue) circle with 
a dashed error bar. Note, as described in Sec. VB, this uncertainty estimate depends on one's choice of Aqcd- In this paper, we 
use Aqcd = 0.7 GeV Choosing Aqcd = 0.5 GeV would cut the error on n in half and decrease the error on riAi. 

Next we estimate the contribution from the dimension-seven operator {icr ■ B, D^}. Using the notation of Ref. [23], summa- 
rized in Sec. A 4, this operator's contribution to the hyperfine spUtting has a coefficient 



1 



1 



(m^/a)'^ (m4a)^' 



(6.4) 



where the equality holds at the tree level for the choices of parameters in our action. The difference between ami and am2 
captures the discretization error. The fractional error in the hyperfine splitting due to this mismatch is 



(aAQCD)^2am2 



1 



1 



(2arn4)3 (2am2)'^ 



(6.5) 



This error is plotted as a (green) dash-dot line on an X in Fig. 8. It would be added in quadrature with the error on the filled dot, 
if it were to be included in the total error. Again we take Aqcd =0.7 GeV, but choosing Aqcd = 0.5 would cut these error bars 
in half. The error from Eq. (6.5) is small for the Dg splitting at the fine lattice spacing, but increasingly large and non-negligible 
at the coarse and medium-coarse lattice spacings; for the Bg splitting, the error is negligible. 

Finally, we turn to the effects of the 0{as) mistuning in cb, which leads to an 0(as) mismatch between msa and m2a. 
Ideally, cb should be adjusted so the coefficient of hl^^Hcr ■ Bh^~^'> equals Zb^vh^, where Zb is a coefficient with an anomalous 
dimension, such that Zsh^'^^ia ■ Bh'^'^^ is scale and scheme independent [49]. In practice, cb is chosen in some approximation, 
in our case the tadpole-improved tree level of perturbation theory. 

Given a value of cb, our simulations produce 



Ml - Ml = Ai = 



4Ao 



2mB{cB) 



(6.6) 



From Eq. (A26), we see that l/ams has a contribution cb /(I + moa). Hence, to include the leading correction to the hyperfine 
splitting, we shift 



4A2a 



1 




4A2a 


_2amB{cB). 



+ 



dt^^ - CB 



2amB{cB) 2{1 + mQa) 



(6.7) 



Other dimension-six and -seven operators are either redimdant, loop-suppressed, or known to have small coefficients [23]. 
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FIG. 7: The hyperflne splitting, in units of ri, versus the ratio of the light to strange sea-quark masses mi/m's on (a) fine and (b) coarse 
ensembles. Errors are the average 68% bootstrap error. The upper panel in each plot is for charm-like splittings and the lower panel is for 
bottom-like splittings. Values of k are 0.127, 0.0923 for the fine ensembles and 0.122, 0.086 for coarse ensembles. Values of am'g are 0.0272 
and 0.0415 for the fine and coarse ensembles, respectively. 



where Cg*^*^' is the ideal choice. (Because loop corrections to l/ams depend on cb, subleading corrections also exist.) To 
estimate the error in Ai, we have to estimate c^'''^' — cb- In fact, Eq. (6.7) can also be used to shift the central value of the 
hyperfine splitting. 

Reference [50] describes preliminary work on a calculation of the one-loop corrections to c^' , as a function of the bare quark 
mass. For all relevant values of ra^a, the one-loop effects are a small correction to the tadpole-improved Ansatz cb = ih)^, 
provided that ?io is the average link in Landau gauge. On the coarse ensembles, we chose mq this way, and we can estimate 
the remaining correction directly from the calculation in Ref. [50]. Given further uncertainties from higher orders, we take this 
small correction as an uncertainty estimate. On the medium-coarse and fine ensembles, however, we chose Uq to be the average 
plaquette. In those cases, the leading correction to cb comes from. 



^ideal 



Cb 



^O.LL 



U, 



'0,plaq 



(6.8) 



where the labels refer to "Landau-gauge link" and "plaquette." Equation (6.8) leads to significant corrections to the hyperfine 
splitting, so we shift Ai on the medium-coarse and fine ensembles by the amount corresponding to Eq. (6.7) and (6.8). These 
shifts put Ai at the medium-coarse and fine lattice spacings on the same footing as those at the coarse spacing. Empirically, they 
flatten the lattice-spacing dependence. 

For the medium-coarse and fine data, we use the values of wq given in Table XVI to calculate the shift described above. It is 
displayed in Fig. 8 as a (pink) star with a single-sided, positive error bar. To obtain an error bar corresponding to the one-loop 
correction to cb in Ref- [50], we take as(0.09 fm) = 1/3 and use one-loop running to obtain values of as for the coarse and 
medium-coarse lattices. These corrections are shown in in Fig. 8 as a (red) triangle with a solid error bar. 

In summary, discretization errors in the hyperfine splitting are small at the fine lattice spacing; therefore, we take as our final 
results the splittings calculated on the fine lattice. In addition, since the effect of the leading 0{as) mistuning of cb can be 
quantified, we shift our final central values by this amount. All other discretization errors are included in our final error. We 
convert our results to physical units using the values of ri /a and r\ as Usted in Table I. After including the error from the units 
conversion in the total, our final results for the hyperfine splittings are 



Ad, = 145 ± 15 MeV 
As = 40 ± 9 MeV 



(6.9) 
(6.10) 



These results are in good agreement with the PDG values of 143.9 ± 0.4 MeV and 46.1 ±1.5 MeV, respectively. 



C. The critical hopping parameter Kcrit 



In principle, it is possible to carry out a suite of nonperturbative heavy-quark calculations without knowing ^crit^ but in 
practice Kcrit is useful. In particular, it enters the construction of improved bilinear and 4-quark operators via mga Eq. (2.5). It 
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FIG. 8: Hyperfine splittings in ri units versus the squared lattice spacing (fm^) for the (a) Ds meson and (b) Bs meson. Filled (black) dots 
with a solid error bar show the splitting with an error from all sources except discretization. Open (blue) circles with a dashed error bar show 
the splitting with an error that also includes discretization error effects in k. (Green) X's with dash-dotted error bars show the estimated size of 
discretization effects from the lattice-continuum mismatch of the dimension-7 operator {icr ■ B, D^} — the errors are barely visible for the 
Bs system. (Pink) stars with a dotted error bar show the 0{as) discretization error from the 1-loop mismatch between m2 and ms. For the 
difference between the 0{as) discretization effects on the coarse lattice versus the fine and medium-coarse lattices, see the text. 



also enters the computation of matching factors such as Zy and Za [18]. Note that these all amount to small corrections, so we 
do not need a very precise determination of Kcrit- Equation (2.5) shows that it does not have to be much better determined than 
Kc and Kb. 

A nonperturbative definition of Kd-a is the value of k such that the mass of a pseudoscalar meson consisting of two Wilson 
quarks (with the clover action) vanishes. The computation of these light-light pseudoscalar meson masses shares code with the 



TABLE XIV: Percent errors in the hyperfine splitting, aAi, of Ds not including discretization effects. 



Uncertainty 


Fine 


Coarse 


Medium-coarse 


Statistical 


2.2 


1.9 


1.9 


K tuning 


8.8,-7.5) 


(4.0,-3.1) 


(4.0,-2.7) 


Valence nis 











Sea-quark masses 


3.6 


5.4 


6.9 


Total 


(10,-9) 


(7,-7) 


(8,-8) 



TABLE XV: Percent errors in the hyperfine splitting, aAi, of Bs not including discretization effects. 



Uncertainty 


Fine 


Coarse 


Medium-coarse 


Statistical 


9.5 


4.0 


5.6 


K tuning 


(12,-11) 


(17,-17) 


(11,-10) 


Valence iris 











Sea-quark masses 


17 


7.8 


2.6 


Total 


(23,-22) 


(19,-19) 


(13,-12) 



TABLE XVI: Tadpole-improvement factors for the estimate of the 0(qs) discretization error shown in Fig. 8. 



ensemble 


^O.plaqucttc 


^0, Landau 


fine (0.0062, 0.031) 


0.878 


0.854 


medium-coarse (0.0097, 0.0484) 


0.860 


0.822 
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TABLE XVII: Values of Kcrit by ensemble, "mq used" gives the origin of the wo value used in the /tcrit determination. Kan values are given 
in two columns. The first /tcrit column contains values which were determined by a fit. The second /tcrit column contains values which 
were estimated from fitted values at the same (approximate) lattice spacing. The last column gives the fit method used in the determination, 
explained in the text. 



/tcrit 



Lattice 


{am'i, am'g) 


uo used 


Iterated fit Direct fit Estimated 


Fine 


(0.0031,0.031) 
(0.0062, 0.031) 
(0.0062, 0.031) 
(0.0124, 0.031) 


Landau-gauge link 
Landau-gauge link 

plaquette 
Landau-gauge link 


0.1372 
0.1391 
0.1372 




0.1372 


Coarse 


(0.005, 0.050) 

(0.007, 0.050) 
(0.010, 0.050) 
(0.020, 0.050) 
(0.030, 0.050) 


Landau-gauge link 

Landau-gauge link 
Landau-gauge link 
Landau-gauge link 
Landau-gauge link 


0.1379 
0.1378 
0.1377 




0.1379 
0.1379 


Medium-coarse (0.0097, 0.0484) 


plaquette 




0.1424 






(0.0194, 0.0484) 


plaquette 




0.1424 






(0.0290, 0.0484) 


plaquette 




0.1423 





work reported here and in Ref. [8], and it is convenient to report the analysis here. The value of /«crit depends on uo via our 
choice of clover coupling, cb = c-e = u^^- In this and other work [3-8], uq has been set sometimes from the average plaquette 
and sometimes from the average link in Landau gauge. The prescription for uq used in each /tcrit determination is given in 
column four of Table XVII. 

The determination of Kcih is carried out on a subset of the available configurations, 50-100 configurations for the fine en- 
sembles and 400-600 for the coarse and medium-coarse. We compute two-point correlators for a range of k that yields meson 
masses of about Mps = 450-900 MeV on the fine ensembles, 650-1100 MeV on the coarse ensembles, and 550-950 MeV 
on the medium-coarse ensembles. It is impractical to push to lower Mpg due to exceptional configurations. Mps is a func- 
tion of the quark mass, which we parametrize as the tree-level, tadpole-improved kinetic or rest mass. In the relevant region, 
mi_2a = moa[l — |moa] + O{{moa)^), so both pertain equally well. The meson masses can be fit to a polynomial ansatz 

a^M|g(/t) = A + Bam2{K, Kcrit) + Ca^ml{K, /tcrit) (6.11) 

(or TOi instead of 7712), where A = when /tcrit is correctly adjusted. 

We use two techniques to determine /tcrit- One method starts with a reasonable value of /tcrit and fits Eq. (6.11) to obtain A, 
B, and C, which depend implicitly on /tcrit- A better trial value of /tcrit is chosen, and the process is iterated until a /tcrit is 
found such that A = 0. We call this the "iterated fit". The second method freezes A to zero, and then B, C, and /tcrit are the fit 
parameters. We call this the "direct fit". On several ensembles the /tcrit values were simply estimated from the other ensembles 
with the same (approximate) lattice spacing, these are labeled as "estimated". 

Table XVII contains our results for /tcrit, indicating the method used. The table does not include error bars for Kcrit, but we 
believe that the results are correct to the number of significant figures shown, even though the range of Afps is high. We carried 
out several tests to verify this accuracy. We compared linear iterated fits [i.e., C = in Eq. (6.1 1)] to the baseline quadratic. We 
also compared direct fits with and without the (continuum) chiral log. These test show that higher order or log contributions do 
not alter our values of /tcrit significantly. We fit comparable data with staggered valence quarks allowing (TOoa)crit 7^ 0, thereby 
testing whether a range of such large Mps skews the results. None of these tests suggests an error larger than a few in the fourth 
digit. Such errors are negUgible compared to those for /tc and Kb — see Tables XII and XIII — when forming moa with Eq. (2.5). 

VII. SUMMARY AND OUTLOOK 

An accurate and precise determination of Ate and /tb is important for all calculations using the Fermilab action [3-8]. In this 
analysis, the error on /tb is dominated by statistics, and the error on /tc receives approximately equal contributions from statistics 
and discretization effects. These errors play a significant role on quantities as diverse as D- and B-meson decay constants [3] 
and the quarkonium hyperfine splitting [8]. Our final results for /tc and /tb are given in Table Xni. 

Another ingredient that is useful for matrix elements [3-5] is the additive renormalization of the bare quark mass or, equiva- 
lently, /tcrit- The improvement and matching of the operators needed to compute these matrix elements depends mildly on /tcrit 
via moa [18]. Our final results for /tcrit are given in Table XVII. 

The key ingredient needed to determine /tc and /«b is a computation of the pseudoscalar and vector heavy-strange meson 
masses. These can be combined to yield the hyperfine splitting for Dg and Bg mesons. Our final results for the hyperfine 
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splittings are given in Eqs. (6.9) and (6.10). Both are in good agreement with the corresponding PDG averages. These results 
bolster confidence in the tuning of Kc and Kh, as well as the choice cb — u^'^. Further tests of these choices come from related 
calculations of the quarkonium spectrum [8] . With detailed attention given to the connection between action parameters and 
mass spUttings, those results are found to be consistent with experiment within the expected uncertainties. 

Improved determinations of Kc, Kb, and Kcrit for the medium-coarse, coarse, and fine ensembles are underway with higher 
statistics, as well as calculations on the new superfine (a w 0.06 fm), and ultrafine (a « 0.045 fm) lattices. The increased 
statistics will also allow us to use higher momentum data and fit to the O(p^) terms in the dispersion relation. Refinements in the 
determination and use ofri/a are allowing for a better understanding of sea-quark effects which will be needed as the statistical 
error on aM2 decreases. We are also investigating the use of twisted boundary conditions [51] which will allow us to obtain data 
points at lower momenta. 

As uncertainties in M2 and Mi decrease, there will be a need for a better understanding of the chiral behavior of these 
masses. One-loop, 0{A/m.Q) chiral perturbation theory results exist for continuum QCD [47]. The extension to staggered 
chiral perturbation theory should be straightforward, and would allow us to extrapolate the light-valence mass to the physical 
up/down quark mass and determine the hyperfine splittings of the and mesons. In this paper, we have included the 
partially quenched expression for the hyperfine splitting in Appendix D, since it is useful in estimating uncertainties from the 
unphysical sea-quark masses. 

In addition, tuned values of Kc, Kb, and Kcrit combined with one-loop (lattice) perturbation theory can yield determinations 
of the pole masses m\ and m2 for both charmed and bottom quarks^, which can be converted to the potential- subtracted, MS, 
and other schemes [6]. Quark masses combined with staggered chiral perturbation theory for the and mesons, can 
yield ab initio calculations of HQET matrix elements [22, 52], which are used to calculate the Cabibbo-Kobayashi-Maskawa 
matrix element |T4b| via inclusive decay measurements. Finally, improved determinations of the oscillating-state energy 
could make determinations of the experimentally accessible masses of the positive parity states, £)*o(2317) and £>si(2460) [53] 
a viable option [54]. 
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Appendix A: Discretization Error in tlie Kinetic Meson Mass 

In this appendix, we present a semi-quantitative estimation of the discretization error in the kinetic mass of heavy-light 
hadrons. We use a formahsm that appUes when both quarks are non-relativistic, even though this approximation is not good for 
the light quark in a heavy-light meson. A posteriori, we examine two ways to re-interpret the resulting formula for a relativistic 
light quark. Both estimates are numerically the same, so we proceed to use the formula in Sec. V B. 

In what follows, the generaUzed masses m\,m2,mi and the coefficient W4, are used to describe the discretization errors. 
Expressions for them when using the Fermilab action are in Refs. [14] or [23] and are given at the end of this appendix for 
convenience. We assume that the light quark (s) has a mass in lattice units <C 1 and makes no significant contribution to 
the discretization error. 

The bound state's kinetic mass can be read off from its kinetic energy (by definition). It will have a kinematic contribution, 
from the constituents' kinetic energy, and a dynamical contribution, from the interaction that binds the constituents. We consider 
each in turn. 



^ The determination of the quark mass from mi requires a non-perturbative calculation of the binding energy as defined by Mi — mi [6]. 
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1. Contributions from constituents' kinetic energy 



The hadron of interest is a heavy-strange meson, a bound state of a heavy quark Q (momentum Q) and a strange antiquark s 
(momentum s). The non-relativistic kinetic energy is 



T = miQ + 



2m2Q 



Smlg 



■ mis + 



2m2s 



(Al) 



The binding energy is communicated to the bound-state kinetic mass via the terms quartic in the momenta and via corrections 
to the potential, given below. In general, the lattice breaks relativistic invariance, so mi ^ TO4, Wi ^ 0. Re-writing the 

kinetic energy in center-of-mass coordinates 



Q 

s 

p 
p 



I2Q 



ni2Q + 

m2Q + m2s 
Q + s, 

m2sQ - m2QS 

m2Q + m2a 



p+p, 



p-p, 



(A2) 

(A3) 
(A4) 
(A5) 



one finds 



T = miQ + mu + 



P^ tP_ _ P^p^ + 2{P-p) 

2(to2Q -I- m2s) 2/i2 4(m2Q + TO2s)^ 



'•2Q 



m 



4Q 



m2s 



Ms 



{in2Q + m2sf 



{wiQUilQ + Wismls) + ■ 



1 



1 1 

+ . 

m2Q m2s 



(A6) 
(A7) 



The only quartic terms shown are those quadratic in P; the omitted terms are not smaller; they just do not contribute to the 
bound state's kinetic energy. The objective is to collect all terms quadratic in P, because their overall coefficient will yield the 
bound state's kinetic mass. 



2. Contribution from tlie interaction: Breit equation 



To obtain the two-particle system's potential energy, one has to work out the scattering amplitude from one-gluon exchange, 
obtaining an expression called the Breit equation [46, 55]. 
In momentum space, for the color-singlet channel 



V{K) = -CFg^D^,{K)MQ{Q + K)u{(,',Q + K)K^Q{Q + K,Q)u{^,Q)MQ{Q) 
xK{s)v{^, s)A:{s, s - K)v{e, s - K)K{s - K), 



(A8) 



where D^j/ is the (lattice) gluon propagator, is the lattice vertex function (for q = Q, s), and Afq is an external-line factor 

needed with the normahzation conditions on spinors employed here [14, 23]. (In continuum field theory, ^f = -^/m/E.) 
To the accuracy needed here, the gluon propagator can be replaced with the continuum propagator. The heavy-quark line is 



Jq = MQ{Q + K)u{i\Q + K)K%{Q + K,Q)u{^,Q)NQ{Q) 



^ - 2iS -{KxQ) 



«(^,0), 



Jq = Mq{Q + K)u{£,',Q + K)Aq{Q + K,Q)u{^,QWq{Q) 
'Q+\K i-ExK 



m2Q 



2msQ 



(A9) 



(AlO) 



and, to the extent that the strange antiquark is non-relativistic, one has a sinnilar expression for the antiquark line = 

Me{s)v{i, s)K{s, s - K)v{e, s - K)Ms{s - K). 
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In Coulomb gauge, 



D44{K) 



K2' 



if2 



(All) 



and the other components vanish. Thus, noting that K4, = i[{Q + KY — Q'^]/2mQ is subleading. 



V{K) = -Cpg^ 



1 



Q s- 



Q KK s 
10 



1 



+ spin-dependent terms. 



(A12) 



Let us discuss each part of the bracket in turn. The leading term yields, after Fourier transforming to position space, the 1/r 
potential. The second yields a contact term proportional to S{r): it is a relativistic correction to the bound state's rest mass, so it 
is of no further interest here. Similarly, the spin-dependent terms do not contribute to the bound state's kinetic energy, so they 
are not written out. The remaining exhibited contributions do contribute to the bound state's kinetic energy, when Q and s are 
eliminated in favor of P and p. 

Next we Fourier transform from it to r using 



1 W 



1 

47rr ' 



r)3 (^2)2 - 2V"»^ 

Following with the substitution of P and p for Q and s this yields 



J (277 



)3 ■ 



V{r,P,p) = - 



CpOts 



1 



2(m2Q + m2sy 



p p 



r 

^ r 2{m2Q + m2s)2 



(A13) 
(A14) 

(A15) 



where the omitted terms do not influence the bound state's kinetic energy. 

Note that K changes p but not P, so r is conjugate to p. To take expectation values, we use the virial theorem 



inVjVir)) 

so the total energy of the bound state, E{P) = {T + V),is 

(p2) /Cpas 



(PiPj) 

M2 



(A16) 



E{P) 



miQ + mis + 

P2 

2(m2Q +m2s) 

p2 



2/i2 

1 - 

(p2 



+ 



2(to2Q + rn2s)^ 2/^2 

PjPj jpiPj) 
{m2Q + m2sY 2/Lt2 

P!'<pi) 



2m2(to2Q+TO2s) {m2Q+m2s)\ r 



1 - M2 









m\s^ 








m\s 



(m2Q m2sy 

The first line of Eq. (A17) shows the binding energy adding to the quarks' rest masses to form the bound state's rest mass, 

(p2) ICpa 



(A17) 



Ml = miQ -h TOis -f- 



2^2 



(A18) 



The second fine shows the same binding energy modifying the kinetic energy. The remaining terms are discretization errors. In 
general they are a bit messy, but they simphfy for the S'-wave states we use to tune k. Then (ftPj) = 1^,^(2?^), whence 



E{P) = Ml 



2Mo 



(A19) 
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where 



M2 = m2Q + m2s 



+ 



3 



M2 



"•2Q , ™L 



4^3(p!) 
3 2/X2 



The last line exhibits the discretization errors, which would vanish if 7714 = 7712,^4 = 0. 
The error can be re- written 



SMo 



3 2^2 



M2 



^2Q ^ rn^ 



^a^,2 [wiQ{m2Qaf + Wis{m2saf 



(A20) 



(A21) 



which is equivalent to Eq. (14) of Ref. [46]. Note that the error ends up being proportional to the internal kinetic energy of the 
bound state, (p^)/2/i2. 



3. Relativistic light degrees of freedom 

For asqtad light quarks, the discretization errors are O(asm^a^) and O(TO^a^). So, for a semi-quantitative estimate of the 
discretization error, it should be safe to assume nriis = m2s = mis = '^s. a^Wis = 0. Equation (A21) is then 



SMo 



1 ip') 
3m2Q 2^2 



H2 



'■2Q 



'4(3 



1 + ^WiQ{m2Qaf 



(All) 



To use this formula we need a value for (p^), and we consider two possibihties. The first is to replace (p^) with A^. The reduced 
mass H2 then cancels, yielding a sensible limit even when rrig 0. The second is to replace the non-relativistic kinetic energy 
(p^) /2/i2 with a relativistic version, namely A. If we take a constituent quark mass iris = 5 A, then this discretization-error 
estimate equals that of the first approach to 0(ms/mQ). 



4. The generalized masses and Wi 

General tree-level expressions for the quark masses and W4, were originally given in Ref. [14] and succinctly recapitulated in 
Ref. [23]. For convenience we give them here with parameters C = 1 = ''s as in our simulations 

moa = —[^- ) , (A23) 



Uq \2k 2Kcrit, 

m\a = ln(l -|- moa) (A24) 

<A25, 



m2a moa{2 + moa) 1 + moa 
' ' +^^, (A26) 



msa moa(2 -|- moa) 1 -|- moa 

' ' (A27) 



Am^a^ [moa(2 + moa)]^ moa(2 + moa) ' 

1 - 8 I 4 + 8(l + moa) ^ 1 ^^^8) 

m\a'^ [moa(2 -|- moa)]^ [moa(2 -|- moa)]^ (l-|-moa)^' 

2 1 

"^4 = 77. 7 + irr., (A29) 

moa(2 -h moa) 4(1 -|- moa) 



These expressions and Eq. (A22) are used to obtain Table IX. 
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Appendix B: Tables of the kinetic mass 

In this appendix, we tabulate values of the pseudoscalar, vector, and spin-averaged kinetic mass, aM2, aM|, and M2, respec- 
tively. Values are given for all combinations of k and am'g on the ensembles used for tuning k. x^/dof and the p vidue, one 
minus the cumulative distribution [45], from the dispersion relation fits are also given. 



_ xVdof (p) 

K aM2 aM^ 0M2 aM2 0M2 

om; = 0.0272 0.090 2.30(17) 2.31(25) 2.31(21) 0.21 (0.81) 0.24(0.79) 

0.0923 2.19(15) 2.22(22) 2.21(19) 0.22(0.80) 0.35 (0.71) 

0.093 2.16(14) 2.19(22) 2.18(18) 0.22 (0.80) 0.37 (0.69) 

0.1256 0.860(19) 0.936(42) 0.917(32) 0.09(0.92) 0.14(0.87) 

0.127 0.819(15) 0.912(39) 0.889(30) 0.36(0.70) 0.00(1.0) 



om; = 0.031 0.0923 2.22(14) 2.22(21) 2.22(18) 0.18(0.84) 0.31 (0.73) 

0.1256 0.871(18) 0.947(38) 0.928(30) 0.14(0.87) 0.16(0.85) 
0.127 0.828(15) 0.918(37) 0.895(29) 0.40(0.67) 0.00(1.0) 



TABLE XVlll: The kinetic meson mass for bottom- and charm-type mesons on the (0.0062, 0.031) fine ensemble from fits to E^{p) — S^(0) 
using |n| < Fits are done to obtain aM2 and 0M2 and the results are then spin averaged. Uncertainties are the average 68% bootstrap 
error, x^/dof with the p value in parentheses is also given. The p value is one minus the x^ cumulative distribution [45]. 









xVdof (p) 


K aM2 




aM 


aM2 aM| 



om;=0.03 0.074 3.78(49) 3.64(54) 3.67(50) 1.12(0.33) 0.17(0.84) 

0.086 2.93(21) 3.03(33) 3.01(29) 0.28 (0.76) 0.21 (0.81) 

0.093 2.50(14) 2.66(24) 2.62(21) 0.05(0.95) 0.11 (0.90) 

0.119 1.263(16) 1.402(43) 1.368(34) 0.84(0.43) 0.20(0.82) 

0.122 1.132(17) 1.270(46) 1.236(37) 0.35(0.70) 0.34(0.71) 

0.124 1.038(16) 1.161(43) 1.130(33) 0.28(0.76) 0.52(0.60) 

am'g = 0.0415 0.074 3.66(35) 3.75(53) 3.73(48) 0.26 (0.77) 0.23 (0.79) 

0.086 2.99(19) 3.09(28) 3.06(25) 0.46 (0.63) 0.48 (0.62) 

0.093 2.57(13) 2.75(25) 2.70(21) 0.14(0.87) 0.31 (0.73) 

0.119 1.292(15) 1.456(41) 1.415(33) 0.88(0.41) 0.38(0.69) 

0.122 1.157(17) 1.310(44) 1.272(36) 0.19(0.83) 0.24(0.79) 

0.124 1.065(15) 1.200(43) 1.166(34) 0.15(0.86) 0.47(0.63) 



TABLE XIX: Same as Table XVIH but for mesons on the (0.007, 0.050) coarse ensemble. 
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xVdof (p) 


K 


aM2 




aM 


aM2 aM| 


am'q = 0.0484 0.070 


4.54(32) 


4.53(45) 


4.53(41) 


0.55 (0.58) 0.78 (0.46) 


0.080 


3.79(19) 


3.77(27) 


3.78(24) 


0.54(0.58) 1.19(0.31) 


0.115 


1.747(25) 


1.825(47) 


1.805(37) 


1.32(0.27) 0.36(0.70) 


0.125 


1.304(12) 


1.415(32) 


1.387(26) 


1.23 (0.29) 0.05 (0.95) 



am'g = 0.0387 0.070 4.47(35) 4.44(50) 4.44(46) 0.47 (0.63) 0.72 (0.49) 

0.080 3.73(21) 3.70(31) 3.71(28) 0.53 (0.59) 1.07 (0.34) 

0.115 1.725(28) 1.804(58) 1.784(47) 1.06(0.43) 0.04(0.96) 

0.125 1.282(13) 1.388(37) 1.361(29) 0.89(0.41) 0.07 (0.93) 



TABLE XX: Same as Table XVIII but for mesons on the (0.0097, 0.0484) medium-coarse ensemble. 
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Appendix C: Tables of Mi = E{0) and the hyperflne splitting 

In this appendix, we tabulate the hyperflne splitting aAi and ri Ai discussed in Sec. VIB. 

1. The hyperflne splitting in lattice units aAi 

In this subsection, we tabulate values of aAi relevant to the discussion in Sec. VIB of the uncertainty in the hyperflne splitting 
due to statistics, k tuning, and the Ught valence mass. 



ensemble oMi aMi aAi 



am;= 0.0272 0.090 (0.0062,0.031) 1.7387(13) 1.7546(19) 0.0158(15) 

0.0923 (0.0031,0.031) 1.6877(21) 1.7054(25) 0.0177(19) 

0.0923 (0.0062,0.031) 1.6870(13) 1.7037(20) 0.0167(16) 

0.0923 (0.0124,0.031) 1.6835(16) 1.7024(19) 0.0188(14) 

0.1256 (0.0062, 0.031) 0.8408(8) 0.8968(16) 0.0561(15) 

0.127 (0.0031 0.031) 0.7944(9) 0.8534(19) 0.0590(19) 

0.127 (0.0062,0.031) 0.7946(7) 0.8544(15) 0.0599(13) 

0.127 (0.0124,0.031) 0.7901(7) 0.8514(11) 0.0613(12) 

amq= 0.031 0.090 (0.0062,0.031) 1.7441(13) 1.7601(17) 0.0159(14) 

0.0923 (0.0062,0.031) 1.6926(12) 1.7093(18) 0.0167(14) 

0.1256 (0.0062, 0.031) 0.8470(8) 0.9030(14) 0.0560(13) 

0.127 (0.0062, 0.031) 0.8009(7) 0.8606(14) 0.0597(12) 



TABLE XXI: Fine-ensemble values of the rest mass Mi = E{0) and hyperflne splitting Ai. am'g = 0.0272 and 0.031. Uncertainties are the 
average 68% bootstrap error 



K ensemble oMi aMl aAi 

om;= 0.0415 0.074 (0.007,0.050) 2.2394(22) 2.2618(25) 0.0224(09) 

0.086 (0.005,0.050) 1.9662(17) 1.9941(27) 0.0279(18) 

0.086 (0.007,0.050) 1.9644(17) 1.9943(21) 0.0299(11) 

0.086 (0.010,0.050) 1.9676(16) 1.9978(21) 0.0301(12) 

0.086 (0.020, 0.050) 1.9584(16) 1.9891(21) 0.0307(14) 



0.122 (0.005,0.050) 1.0529(10) 1.1399(22) 0.0870(17) 

0.122 (0.007,0.050) 1.0520(7) 1.1393(17) 0.0873(15) 

0.122 (0.010,0.050) 1.0549(10) 1.1414(26) 0.0865(22) 

0.122 (0.020,0.050) 1.0446(9) 1.1339(16) 0.0894(16) 

0.124 (0.007, 0.050) 0.9871(7) 1.0819(17) 0.0948(15) 

om;= 0.030 0.074 (0.007,0.050) 2.2241(26) 2.2466(29) 0.0225(11) 

0.086 (0.007, 0.050) 1.9488(21) 1.9787(25) 0.0299(14) 

0.122 (0.007,0.050) 1.0339(08) 1.1220(20) 0.0881(17) 



TABLE XXII: Same as Table XXI but for the coarse ensembles with am, = 0.0415 and 0.03. 



K ensemble aMi aM^ aAi 

am;= 0.0484 0.076 (0.0097.0.0484) 2.3192(27) 2.3472(36) 0.0280(17) 

0.076 (0.0194,0.0484) 2.3153(30) 2.3424(47) 0.0270(26) 

0.076 (0.0290,0.0484) 2.3137(23) 2.3445(21) 0.0308(16) 

0.080 (0.0097, 0.0484) 2.2298(24) 2.2606(34) 0.0308(17) 

0.122 (0.0097, 0.0484) 1.2427(8) 1.3390(21) 0.0963(19) 

0.122 (0.0194, 0.0484) 1.2397(8) 1.3400(17) 0.1004(14) 

0.122 (0.0290, 0.0484) 1.2364(7) 1.3402(17) 0.1038(14) 

0.125 (0.0097,0.0484) 1.1565(8) 1.2634(21) 0.1069(20) 

om;= 0.0387 0.076 (0.0097,0.0484) 2.3060(31) 2.3341(40) 0.0281(19) 

0.122 (0.0097, 0.0484) 1.2271(9) 1.3237(25) 0.0966(24) 



TABLE XXm: Same as Table XXI but for medium-coarse ensembles with am(, = 0.0484 and 0.0387. 
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2. The hyperflne splitting in physical units ri Ai 



In this subsection, we tabulate values of riAi relevant to the discussion in Sec. VIB of the dependence of the hyperflne 
splitting on the sea-quark masses. 



ensemble riAi 
K = 0.0923 (0.0031,0.031) 0.0653(69) 
(0.0062,0.031) 0.0618(58) 
(0.0124,0.031) 0.0700(52) 

K = 0.127 (0.0031 0.031) 0.2178(69) 
(0.0062, 0.031) 0.2217(48) 
(0.0124, 0.031) 0.2281(43) 



TABLE XXIV: Fine-ensemble values of the hyperflne splitting Ai in units of ri. am'q = 0.0272. Uncertainties are the average 68% bootstrap 
error. 



ensemble riAi 
K = 0.086 (0.005, 0.050) 0.0738(46) 
(0.007, 0.050) 0.0788(29) 
(0.010, 0.050) 0.0788(32) 
(0.020, 0.050) 0.0814(36) 

K = 0.122 (0.005,0.050) 0.2301(44) 
(0.007, 0.050) 0.2300(39) 
(0.010, 0.050) 0.2265(58) 

(0.020. 0.050) 0.2370(43) 



TABLE XXV: Same as Table XXIV but for the coarse-ensembles with am'q = 0.0415. 



ensemble riAi 
K = 0.076 (0.0097, 0.0484) 0.0616(37) 
(0.0194, 0.0484) 0.0603(57) 
(0.0290, 0.0484) 0.0699(37) 

/t = 0.122 (0.0097,0.0484) 0.2117(41) 
(0.0194, 0.0484) 0.2244(30) 
(0.0290, 0.0484) 0.2357(39) 



TABLE XXVI: Same as Table XXIV but for medium-coarse-ensembles with am' = 0.0484. 
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Appendix D: Partially quenched chiral perturbation theory for the heavy-light hyperfine spUtting 

For full (unquenched) QCD, Jenkins [47] has calculated the hyperfine spUtting at one loop in heavy-meson chiral perturbation 
theory. It is not difficult to take her result (Eq. (A. 10) of Ref [47]) and extend it to partially quenched QCD. The further step 
of including staggered taste- violations (i.e., doing staggered chiral perturbation theory) would also be fairly straightforward, but 
we do not take it here because the continuum partially quenched form is sufficient for estimating the small systematic effect due 
to the mistuning of sea quark masses. Unlike Jenkins, we neglect electromagnetic and isospin-violating effects. 

At the quark-flow level, the relevant diagrams are the self-energy diagrams shown in Fig. 5(a) [left] of Ref. [56] (the "con- 
nected diagram") and in Figs. 5(b), (c) [left] of Ref. [56] (the "disconnected diagram").* One simply needs to determine how 
much of Jenkins's result comes from each of these two diagrams. This is accomplished by noting that, when the Ught valence 
quark is a m (a = 1 in Jenkins's notation), an internal kaon only appears in the connected diagram, when the quark in the virtual 
loop is an s. This fixes the normalization of the connected diagram. Using the methods described in Refs. 1 56-58] (but dropping 
the taste violations and indeed the taste degree of freedom itself), the disconnected diagram is easily calculated. Its normalization 
can then be fixed so that it supplies the remainder of the a = 1 result in Ref. [47]. 

There are ample checks of this reasoning. First, the same normalizations must apply for any choice of the valence mass. The 
r] contributes in each case only through the disconnected diagram, while the pion contributions come from both connected and 
disconnected diagrams for valence u or d (a = 1,2), and must be absent for valence s (a = 3). Finally the contribution from the 
unphysical ss state, which appears in each diagram for a = 3, should cancel. 

It is then immediate to write down the partially quenched version. Let the light valence quark be x, with mass nix, and let 
the sea quarks be u, d, s with masses m„ = nid = mi and nis- With the light meson decay constant / normalized so that 
f K. K. 130 MeV, the hyperfine splitting M* - is given by 

M* - = A - ^^'^log + 2A('^)(2mi + m,) + 2A(°)m^ , (Dl) 

where A is the splitting in the (three-flavor) chiral limit, and A^'^^ and A^"' are LECs that start at order 1 /mg in the heavy quark 
expansion. The non-analytic chiral logarithms ^log are 

5iog= Yl m!F)-lR^x''\{rn}Af^})iiMh-lY.Df;^\{rn},{f^})i{M^). (D2) 

F=u,d,s j=X,ri 

Here Mx is the mass of the valence xx meson, and M^f is the mass of the mixed valence-sea xF meson. The residue functions 
i?^"'*^^ and Dj^-'^\ as well as the chiral logarithm functions (.{m^) and i{vn?), are defined in Refs. [57, 58]. The term with the 
sum over F comes from the cormected diagram, while those with the residue functions come from the discormected diagram, 
which has a double pole at in the partially quenched case. The denominator ({m}) and numerator ({/u}) mass sets are 

{m} = {Mx, , {/x} = {Mu, Ms} (D3) 

with Mu and Ms the masses of the uu and ss mesons, respectively. 
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